WGCNA Analysis: Inv4m DEG Networks

Overview

Apply WGCNA to Inv4m DEGs to identify co-expression modules.

Approach: Unsigned networks with power = 12, minModuleSize = 10, better suited for identifying regulatory hubs like JMJ genes that can both activate and repress targets.

Libraries

library(SummarizedExperiment)
library(tidyverse)
library(WGCNA)
library(pheatmap)

Load Data

Normalized Expression

voomR <- readRDS(file.path(paths$intermediate, "normalized_expression_voom_object_leaf_trt.rds"))

cat("Expression data loaded:\n")

## Expression data loaded:

cat("  Genes:", nrow(voomR$E), "\n")

##   Genes: 24249

cat("  Samples:", ncol(voomR$E), "\n")

##   Samples: 43

DEG Results

effects_df <- read_csv(
  file.path(paths$intermediate, "predictor_effects_leaf_interaction_model.csv")
)

inv4m_degs <- effects_df %>%
  filter(
    predictor == "Inv4m",
    is_DEG == TRUE
  ) %>%
  pull(gene) %>%
  unique()

cat("Inv4m DEGs identified:", length(inv4m_degs), "\n")

## Inv4m DEGs identified: 776

cat("  Upregulated:", 
    sum(effects_df$predictor == "Inv4m" & 
        effects_df$is_DEG & 
        effects_df$regulation == "Upregulated"), "\n")

##   Upregulated: 237

cat("  Downregulated:", 
    sum(effects_df$predictor == "Inv4m" & 
        effects_df$is_DEG & 
        effects_df$regulation == "Downregulated"), "\n")

##   Downregulated: 539

Prepare Data

Filter Expression Matrix

expr_matrix_degs <- t(voomR$E[inv4m_degs, ])

cat("Expression matrix dimensions:", dim(expr_matrix_degs), "\n")

## Expression matrix dimensions: 43 776

cat("  Samples (rows):", nrow(expr_matrix_degs), "\n")

##   Samples (rows): 43

cat("  Genes (columns):", ncol(expr_matrix_degs), "\n")

##   Genes (columns): 776

Sample Metadata

sample_table <- voomR$targets %>%
  as.data.frame() %>%
  rownames_to_column("Sample") %>%
  select(
    Sample,
    Genotype,
    leaf_tissue,
    Treatment
  ) %>%
  mutate(
    Genotype = factor(Genotype, levels = c("CTRL", "Inv4m")),
    leaf_tissue = factor(leaf_tissue)
  )

rownames(sample_table) <- sample_table$Sample

cat("Sample distribution:\n")

## Sample distribution:

with(sample_table, {
  cat("Genotype × Leaf tissue:\n")
  print(table(Genotype, leaf_tissue))
  cat("\nGenotype × Treatment:\n")
  print(table(Genotype, Treatment))
})

## Genotype × Leaf tissue:
##         leaf_tissue
## Genotype 1 2 3 4
##    CTRL  4 6 5 5
##    Inv4m 7 5 6 5
## 
## Genotype × Treatment:
##         Treatment
## Genotype +P -P
##    CTRL  12  8
##    Inv4m 12 11

Network Construction

Combined Network (All Samples)

datExpr <- expr_matrix_degs

cat("Expression data prepared for WGCNA:\n")

## Expression data prepared for WGCNA:

cat("  Samples (rows):", nrow(datExpr), "\n")

##   Samples (rows): 43

cat("  Genes (columns):", ncol(datExpr), "\n")

##   Genes (columns): 776

powers <- c(1:20)

sft_all <- pickSoftThreshold(
  datExpr,
  powerVector = powers,
  networkType = "unsigned",
  verbose = 5
)

## pickSoftThreshold: will use block size 776.
##  pickSoftThreshold: calculating connectivity for given powers...
##    ..working on genes 1 through 776 of 776

##    Power SFT.R.sq   slope truncated.R.sq mean.k. median.k. max.k.
## 1      1  0.49400  2.6700          0.936  324.00   314.000  468.0
## 2      2  0.23300  0.8980          0.813  172.00   161.000  309.0
## 3      3  0.00295 -0.0599          0.672  102.00    93.000  222.0
## 4      4  0.21500 -0.4490          0.671   66.10    58.000  168.0
## 5      5  0.58600 -0.7990          0.716   45.20    37.500  134.0
## 6      6  0.76900 -1.0200          0.792   32.40    25.300  112.0
## 7      7  0.86000 -1.1400          0.855   24.10    17.600   96.4
## 8      8  0.91700 -1.2000          0.902   18.40    12.600   84.5
## 9      9  0.93200 -1.2400          0.915   14.50     8.920   75.2
## 10    10  0.93200 -1.2500          0.912   11.60     6.490   67.7
## 11    11  0.95700 -1.2700          0.946    9.53     4.810   61.5
## 12    12  0.93200 -1.2700          0.915    7.91     3.580   56.3
## 13    13  0.94400 -1.2600          0.935    6.66     2.690   51.8
## 14    14  0.94100 -1.2600          0.929    5.67     2.070   47.9
## 15    15  0.96600 -1.2400          0.961    4.88     1.560   44.5
## 16    16  0.96700 -1.2200          0.964    4.23     1.230   41.5
## 17    17  0.95600 -1.2300          0.950    3.70     0.936   38.7
## 18    18  0.96200 -1.2100          0.961    3.26     0.739   36.3
## 19    19  0.95700 -1.2100          0.951    2.89     0.581   34.1
## 20    20  0.96700 -1.1900          0.963    2.57     0.465   32.1

par(mfrow = c(1, 2))
plot(
  sft_all$fitIndices[, 1],
  -sign(sft_all$fitIndices[, 3]) * sft_all$fitIndices[, 2],
  xlab = "Soft Threshold (power)",
  ylab = "Scale Free Topology Model Fit, signed R^2",
  main = "Scale Independence - All Samples",
  type = "n"
)
text(
  sft_all$fitIndices[, 1],
  -sign(sft_all$fitIndices[, 3]) * sft_all$fitIndices[, 2],
  labels = powers,
  cex = 0.9,
  col = "red"
)
abline(h = 0.80, col = "blue")

plot(
  sft_all$fitIndices[, 1],
  sft_all$fitIndices[, 5],
  xlab = "Soft Threshold (power)",
  ylab = "Mean Connectivity",
  main = "Mean Connectivity - All Samples",
  type = "n"
)
text(
  sft_all$fitIndices[, 1],
  sft_all$fitIndices[, 5],
  labels = powers,
  cex = 0.9,
  col = "red"
)

pow_all <- sft_all$powerEstimate
if (is.na(pow_all)) {
  pow_all <- sft_all$fitIndices %>%
    filter(-sign(Slope) * SFT.R.sq > 0.80) %>%
    slice_head(n = 1) %>%
    pull(Power)
  
  if (length(pow_all) == 0) {
    pow_all <- 12
    cat("\nWarning: No power achieved R^2 > 0.80, using default power = 12\n")
  }
}

cat("\nSelected soft-thresholding power (all samples):", pow_all, "\n")

## 
## Selected soft-thresholding power (all samples): 7

cat("Scale-free topology fit R^2:", 
    sft_all$fitIndices[sft_all$fitIndices$Power == pow_all, "SFT.R.sq"], "\n")

## Scale-free topology fit R^2: 0.8601684

cat("Building UNSIGNED network with blockwiseModules:\n")

## Building UNSIGNED network with blockwiseModules:

net_all <- blockwiseModules(
  datExpr,
  power = pow_all,
  networkType = "unsigned",
  TOMType = "unsigned",
  minModuleSize = 10,
  maxBlockSize = 25000,
  reassignThreshold = 0,
  mergeCutHeight = 0,
  numericLabels = TRUE,
  pamRespectsDendro = FALSE,
  deepSplit = 4,
  verbose = 3
)

##  Calculating module eigengenes block-wise from all genes
##    Flagging genes and samples with too many missing values...
##     ..step 1
##  ..Working on block 1 .
##     TOM calculation: adjacency..
##     ..will not use multithreading.
##      Fraction of slow calculations: 0.000000
##     ..connectivity..
##     ..matrix multiplication (system BLAS)..
##     ..normalization..
##     ..done.
##  ....clustering..
##  ....detecting modules..
##  ....calculating module eigengenes..
##  ....checking kME in modules..
##  ..merging modules that are too close..
##      mergeCloseModules: Merging modules whose distance is less than 0
##        Calculating new MEs...

wgcna_all <- list(
  combined = structure(
    list(
      datExpr = datExpr,
      moduleColors = net_all$colors,
      MEs = net_all$MEs,
      sampleTable = sample_table
    ),
    class = "WGCNA"
  )
)

cat("\nNetwork built successfully\n")

## 
## Network built successfully

cat("  Modules:", length(unique(net_all$colors)), "\n")

##   Modules: 14

cat("  Genes:", length(net_all$colors), "\n")

##   Genes: 776

Genotype-Specific Networks

ctrl_samples <- rownames(sample_table)[sample_table$Genotype == "CTRL"]
inv4m_samples <- rownames(sample_table)[sample_table$Genotype == "Inv4m"]

datExpr_ctrl <- datExpr[ctrl_samples, ]
datExpr_inv4m <- datExpr[inv4m_samples, ]

cat("CTRL expression data:\n")

## CTRL expression data:

cat("  Samples:", nrow(datExpr_ctrl), "\n")

##   Samples: 20

cat("  Genes:", ncol(datExpr_ctrl), "\n")

##   Genes: 776

cat("\nInv4m expression data:\n")

## 
## Inv4m expression data:

cat("  Samples:", nrow(datExpr_inv4m), "\n")

##   Samples: 23

cat("  Genes:", ncol(datExpr_inv4m), "\n")

##   Genes: 776

sft_ctrl <- pickSoftThreshold(
  datExpr_ctrl,
  powerVector = powers,
  networkType = "unsigned",
  verbose = 5
)

## pickSoftThreshold: will use block size 776.
##  pickSoftThreshold: calculating connectivity for given powers...
##    ..working on genes 1 through 776 of 776
##    Power SFT.R.sq  slope truncated.R.sq mean.k. median.k. max.k.
## 1      1   0.4840  1.140          0.944  254.00   256.000  373.0
## 2      2   0.0303 -0.125          0.832  121.00   116.000  235.0
## 3      3   0.4830 -0.555          0.849   68.40    60.900  166.0
## 4      4   0.7670 -0.764          0.930   42.90    35.500  124.0
## 5      5   0.8230 -0.874          0.927   28.80    21.700   96.4
## 6      6   0.8610 -0.949          0.929   20.40    13.800   76.9
## 7      7   0.9070 -0.985          0.953   14.90     9.020   62.5
## 8      8   0.9250 -1.020          0.962   11.30     6.320   51.6
## 9      9   0.9500 -1.060          0.968    8.72     4.430   43.3
## 10    10   0.9400 -1.070          0.960    6.87     3.170   37.0
## 11    11   0.9420 -1.110          0.962    5.50     2.340   32.0
## 12    12   0.9510 -1.120          0.976    4.47     1.720   27.9
## 13    13   0.9450 -1.150          0.971    3.67     1.310   24.5
## 14    14   0.9520 -1.160          0.982    3.05     1.010   21.7
## 15    15   0.9540 -1.180          0.979    2.56     0.781   19.2
## 16    16   0.9250 -1.210          0.956    2.16     0.612   17.1
## 17    17   0.9320 -1.220          0.964    1.84     0.474   15.3
## 18    18   0.9420 -1.220          0.977    1.58     0.369   13.8
## 19    19   0.9450 -1.230          0.977    1.37     0.294   12.4
## 20    20   0.9290 -1.240          0.948    1.19     0.232   11.2

sft_inv4m <- pickSoftThreshold(
  datExpr_inv4m,
  powerVector = powers,
  networkType = "unsigned",
  verbose = 5
)

## pickSoftThreshold: will use block size 776.
##  pickSoftThreshold: calculating connectivity for given powers...
##    ..working on genes 1 through 776 of 776
##    Power SFT.R.sq  slope truncated.R.sq mean.k. median.k. max.k.
## 1      1   0.5030  2.060          0.881 225.000  230.0000 310.00
## 2      2   0.0618  0.329          0.749  95.400   97.1000 167.00
## 3      3   0.0678 -0.277          0.749  48.500   47.7000 103.00
## 4      4   0.3080 -0.681          0.782  27.500   25.9000  68.60
## 5      5   0.4650 -0.870          0.842  16.800   15.0000  48.20
## 6      6   0.5940 -1.000          0.896  10.900    9.1800  35.20
## 7      7   0.6710 -1.080          0.932   7.370    5.8300  26.60
## 8      8   0.7240 -1.170          0.938   5.160    3.8400  20.50
## 9      9   0.7650 -1.190          0.947   3.720    2.6000  16.20
## 10    10   0.7730 -1.250          0.931   2.740    1.7800  12.90
## 11    11   0.8030 -1.260          0.943   2.070    1.2500  10.50
## 12    12   0.8040 -1.300          0.925   1.590    0.8930   8.66
## 13    13   0.8290 -1.310          0.938   1.240    0.6550   7.21
## 14    14   0.8480 -1.310          0.935   0.981    0.4800   6.06
## 15    15   0.8860 -1.280          0.942   0.786    0.3660   5.14
## 16    16   0.8240 -1.330          0.855   0.638    0.2680   4.48
## 17    17   0.8920 -1.410          0.955   0.523    0.2000   4.34
## 18    18   0.2720 -2.640          0.208   0.434    0.1500   4.22
## 19    19   0.3050 -3.650          0.257   0.363    0.1180   4.12
## 20    20   0.3200 -3.670          0.275   0.306    0.0889   4.02

par(mfrow = c(2, 2))

plot(
  sft_ctrl$fitIndices[, 1],
  -sign(sft_ctrl$fitIndices[, 3]) * sft_ctrl$fitIndices[, 2],
  xlab = "Soft Threshold (power)",
  ylab = "Scale Free Topology Fit R^2",
  main = "Scale Independence - CTRL",
  type = "n"
)
text(
  sft_ctrl$fitIndices[, 1],
  -sign(sft_ctrl$fitIndices[, 3]) * sft_ctrl$fitIndices[, 2],
  labels = powers,
  cex = 0.9,
  col = "red"
)
abline(h = 0.80, col = "blue")

plot(
  sft_ctrl$fitIndices[, 1],
  sft_ctrl$fitIndices[, 5],
  xlab = "Soft Threshold (power)",
  ylab = "Mean Connectivity",
  main = "Mean Connectivity - CTRL",
  type = "n"
)
text(
  sft_ctrl$fitIndices[, 1],
  sft_ctrl$fitIndices[, 5],
  labels = powers,
  cex = 0.9,
  col = "red"
)

plot(
  sft_inv4m$fitIndices[, 1],
  -sign(sft_inv4m$fitIndices[, 3]) * sft_inv4m$fitIndices[, 2],
  xlab = "Soft Threshold (power)",
  ylab = "Scale Free Topology Fit R^2",
  main = "Scale Independence - Inv4m",
  type = "n"
)
text(
  sft_inv4m$fitIndices[, 1],
  -sign(sft_inv4m$fitIndices[, 3]) * sft_inv4m$fitIndices[, 2],
  labels = powers,
  cex = 0.9,
  col = "red"
)
abline(h = 0.80, col = "blue")

plot(
  sft_inv4m$fitIndices[, 1],
  sft_inv4m$fitIndices[, 5],
  xlab = "Soft Threshold (power)",
  ylab = "Mean Connectivity",
  main = "Mean Connectivity - Inv4m",
  type = "n"
)
text(
  sft_inv4m$fitIndices[, 1],
  sft_inv4m$fitIndices[, 5],
  labels = powers,
  cex = 0.9,
  col = "red"
)

pow_ctrl <- sft_ctrl$powerEstimate
if (is.na(pow_ctrl)) {
  pow_ctrl <- sft_ctrl$fitIndices %>%
    filter(-sign(Slope) * SFT.R.sq > 0.80) %>%
    slice_head(n = 1) %>%
    pull(Power)
  
  if (length(pow_ctrl) == 0) {
    pow_ctrl <- 12
    cat("\nWarning: CTRL - No power achieved R^2 > 0.80, using default = 12\n")
  }
}

pow_inv4m <- sft_inv4m$powerEstimate
if (is.na(pow_inv4m)) {
  pow_inv4m <- sft_inv4m$fitIndices %>%
    filter(-sign(Slope) * SFT.R.sq > 0.80) %>%
    slice_head(n = 1) %>%
    pull(Power)
  
  if (length(pow_inv4m) == 0) {
    pow_inv4m <- 12
    cat("\nWarning: Inv4m - No power achieved R^2 > 0.80, using default = 12\n")
  }
}

cat("\nSelected powers:\n")

## 
## Selected powers:

cat("  CTRL:", pow_ctrl, 
    "(R^2 =", sft_ctrl$fitIndices[sft_ctrl$fitIndices$Power == pow_ctrl, "SFT.R.sq"], ")\n")

##   CTRL: 6 (R^2 = 0.8605398 )

cat("  Inv4m:", pow_inv4m, 
    "(R^2 =", sft_inv4m$fitIndices[sft_inv4m$fitIndices$Power == pow_inv4m, "SFT.R.sq"], ")\n")

##   Inv4m: 15 (R^2 = 0.8863985 )

cat("Building CTRL network:\n")

## Building CTRL network:

net_ctrl <- blockwiseModules(
  datExpr_ctrl,
  power = pow_ctrl,
  networkType = "unsigned",
  TOMType = "unsigned",
  minModuleSize = 10,
  maxBlockSize = 25000,
  reassignThreshold = 0,
  mergeCutHeight = 0,
  numericLabels = TRUE,
  pamRespectsDendro = FALSE,
  deepSplit = 4,
  verbose = 3
)

##  Calculating module eigengenes block-wise from all genes
##    Flagging genes and samples with too many missing values...
##     ..step 1
##  ..Working on block 1 .
##     TOM calculation: adjacency..
##     ..will not use multithreading.
##      Fraction of slow calculations: 0.000000
##     ..connectivity..
##     ..matrix multiplication (system BLAS)..
##     ..normalization..
##     ..done.
##  ....clustering..
##  ....detecting modules..
##  ....calculating module eigengenes..
##  ....checking kME in modules..
##      ..removing 6 genes from module 1 because their KME is too low.
##      ..removing 11 genes from module 2 because their KME is too low.
##      ..removing 1 genes from module 3 because their KME is too low.
##      ..removing 1 genes from module 4 because their KME is too low.
##      ..removing 3 genes from module 8 because their KME is too low.
##      ..removing 3 genes from module 10 because their KME is too low.
##      ..removing 1 genes from module 15 because their KME is too low.
##  ..merging modules that are too close..
##      mergeCloseModules: Merging modules whose distance is less than 0
##        Calculating new MEs...

wgcna_ctrl <- list(
  ctrl = structure(
    list(
      datExpr = datExpr_ctrl,
      moduleColors = net_ctrl$colors,
      MEs = net_ctrl$MEs,
      sampleTable = sample_table[ctrl_samples, ]
    ),
    class = "WGCNA"
  )
)

cat("\nCTRL network built\n")

## 
## CTRL network built

cat("  Modules:", length(unique(net_ctrl$colors)), "\n")

##   Modules: 19

cat("  Genes:", length(net_ctrl$colors), "\n")

##   Genes: 776

cat("Building Inv4m network:\n")

## Building Inv4m network:

net_inv4m <- blockwiseModules(
  datExpr_inv4m,
  power = pow_inv4m,
  networkType = "unsigned",
  TOMType = "unsigned",
  minModuleSize = 10,
  maxBlockSize = 25000,
  reassignThreshold = 0,
  mergeCutHeight = 0,
  numericLabels = TRUE,
  pamRespectsDendro = FALSE,
  deepSplit = 4,
  verbose = 3
)

##  Calculating module eigengenes block-wise from all genes
##    Flagging genes and samples with too many missing values...
##     ..step 1
##  ..Working on block 1 .
##     TOM calculation: adjacency..
##     ..will not use multithreading.
##      Fraction of slow calculations: 0.000000
##     ..connectivity..
##     ..matrix multiplication (system BLAS)..
##     ..normalization..
##     ..done.
##  ....clustering..
##  ....detecting modules..
##  ....calculating module eigengenes..
##  ....checking kME in modules..
##  ..merging modules that are too close..
##      mergeCloseModules: Merging modules whose distance is less than 0
##        Calculating new MEs...

wgcna_inv4m <- list(
  inv4m = structure(
    list(
      datExpr = datExpr_inv4m,
      moduleColors = net_inv4m$colors,
      MEs = net_inv4m$MEs,
      sampleTable = sample_table[inv4m_samples, ]
    ),
    class = "WGCNA"
  )
)

cat("\nInv4m network built\n")

## 
## Inv4m network built

cat("  Modules:", length(unique(net_inv4m$colors)), "\n")

##   Modules: 19

cat("  Genes:", length(net_inv4m$colors), "\n")

##   Genes: 776

Network Summary

n_modules_all <- length(unique(net_all$colors))
module_sizes_all <- table(net_all$colors)

cat("\nCombined network (all samples):\n")

## 
## Combined network (all samples):

cat("  Total modules:", n_modules_all, "\n")

##   Total modules: 14

cat("  Total genes:", length(net_all$colors), "\n")

##   Total genes: 776

cat("\nModule sizes:\n")

## 
## Module sizes:

print(sort(module_sizes_all, decreasing = TRUE))

## 
##   1   2   3   4   5   6   7   8   9  10  11  12  13  14 
## 340  77  60  50  35  31  31  30  29  29  21  20  13  10

saveRDS(wgcna_all, file.path(paths$intermediate, "wgcna_all.rds"))
cat("\nNetwork saved to wgcna_all.rds\n")

## 
## Network saved to wgcna_all.rds

n_modules_ctrl <- length(unique(net_ctrl$colors))
module_sizes_ctrl <- table(net_ctrl$colors)

n_modules_inv4m <- length(unique(net_inv4m$colors))
module_sizes_inv4m <- table(net_inv4m$colors)

cat("\nCTRL network:\n")

## 
## CTRL network:

cat("  Total modules:", n_modules_ctrl, "\n")

##   Total modules: 19

cat("  Total genes:", length(net_ctrl$colors), "\n")

##   Total genes: 776

cat("\nModule sizes:\n")

## 
## Module sizes:

print(sort(module_sizes_ctrl, decreasing = TRUE))

## 
##   1   2   3   4   5   6   7   8   9  10  11   0  12  13  14  15  16  17  18 
## 116 102  82  63  56  40  39  36  36  30  28  27  27  18  17  16  15  15  13

cat("\n\nInv4m network:\n")

## 
## 
## Inv4m network:

cat("  Total modules:", n_modules_inv4m, "\n")

##   Total modules: 19

cat("  Total genes:", length(net_inv4m$colors), "\n")

##   Total genes: 776

cat("\nModule sizes:\n")

## 
## Module sizes:

print(sort(module_sizes_inv4m, decreasing = TRUE))

## 
##   0   1   2   3   4   5   6   7   8   9  10  11  12  13  14  15  16  17  18 
## 288  62  50  45  44  40  38  33  26  21  19  17  16  15  14  13  12  12  11

saveRDS(wgcna_ctrl, file.path(paths$intermediate, "wgcna_ctrl.rds"))
saveRDS(wgcna_inv4m, file.path(paths$intermediate, "wgcna_inv4m.rds"))
cat("\nGenotype-specific networks saved\n")

## 
## Genotype-specific networks saved

Differential Module Expression

MEs_all <- net_all$MEs

ME_data_all <- MEs_all %>%
  as.data.frame() %>%
  rownames_to_column("Sample") %>%
  left_join(
    sample_table %>% select(Sample, Genotype, leaf_tissue, Treatment),
    by = "Sample"
  )

module_names_all <- grep("^ME", colnames(MEs_all), value = TRUE)

dme_results <- map_dfr(module_names_all, function(mod) {
  me_values <- ME_data_all[[mod]]
  genotype <- ME_data_all$Genotype
  
  inv4m_vals <- me_values[genotype == "Inv4m"]
  ctrl_vals <- me_values[genotype == "CTRL"]
  
  t_result <- t.test(inv4m_vals, ctrl_vals)
  
  tibble(
    module = mod,
    mean_Inv4m = mean(inv4m_vals),
    mean_CTRL = mean(ctrl_vals),
    diff = mean_Inv4m - mean_CTRL,
    t_statistic = t_result$statistic,
    p_value = t_result$p.value
  )
})

dme_results <- dme_results %>%
  mutate(fdr = p.adjust(p_value, method = "fdr")) %>%
  arrange(p_value)

cat("Differential module expression results:\n")

## Differential module expression results:

print(dme_results)

## # A tibble: 14 × 7
##    module mean_Inv4m mean_CTRL   diff t_statistic  p_value      fdr
##    <chr>       <dbl>     <dbl>  <dbl>       <dbl>    <dbl>    <dbl>
##  1 ME1       -0.141     0.162  -0.303      -51.5  7.30e-39 1.02e-37
##  2 ME5       -0.100     0.115  -0.215       -6.33 1.61e- 7 1.13e- 6
##  3 ME9       -0.0988    0.114  -0.212       -6.08 5.76e- 7 2.69e- 6
##  4 ME7       -0.0909    0.104  -0.195       -5.28 5.39e- 6 1.70e- 5
##  5 ME6        0.0935   -0.107   0.201        5.42 6.08e- 6 1.70e- 5
##  6 ME13      -0.0838    0.0964 -0.180       -4.69 3.11e- 5 7.27e- 5
##  7 ME14      -0.0834    0.0959 -0.179       -4.60 4.32e- 5 8.65e- 5
##  8 ME10      -0.0821    0.0944 -0.177       -4.45 8.25e- 5 1.44e- 4
##  9 ME8        0.0778   -0.0894  0.167        4.24 1.25e- 4 1.95e- 4
## 10 ME3       -0.0774    0.0890 -0.166       -4.04 3.16e- 4 4.42e- 4
## 11 ME11      -0.0737    0.0848 -0.159       -3.77 6.82e- 4 8.68e- 4
## 12 ME4       -0.0606    0.0697 -0.130       -2.92 6.49e- 3 7.19e- 3
## 13 ME2       -0.0597    0.0686 -0.128       -2.89 6.68e- 3 7.19e- 3
## 14 ME12       0.0563   -0.0648  0.121        2.71 1.04e- 2 1.04e- 2

sig_modules <- dme_results %>%
  filter(fdr < 0.05)

cat("\nSignificant genotype-associated modules (FDR < 0.05):", 
    nrow(sig_modules), "\n")

## 
## Significant genotype-associated modules (FDR < 0.05): 14

if (nrow(sig_modules) > 0) {
  print(sig_modules)
}

## # A tibble: 14 × 7
##    module mean_Inv4m mean_CTRL   diff t_statistic  p_value      fdr
##    <chr>       <dbl>     <dbl>  <dbl>       <dbl>    <dbl>    <dbl>
##  1 ME1       -0.141     0.162  -0.303      -51.5  7.30e-39 1.02e-37
##  2 ME5       -0.100     0.115  -0.215       -6.33 1.61e- 7 1.13e- 6
##  3 ME9       -0.0988    0.114  -0.212       -6.08 5.76e- 7 2.69e- 6
##  4 ME7       -0.0909    0.104  -0.195       -5.28 5.39e- 6 1.70e- 5
##  5 ME6        0.0935   -0.107   0.201        5.42 6.08e- 6 1.70e- 5
##  6 ME13      -0.0838    0.0964 -0.180       -4.69 3.11e- 5 7.27e- 5
##  7 ME14      -0.0834    0.0959 -0.179       -4.60 4.32e- 5 8.65e- 5
##  8 ME10      -0.0821    0.0944 -0.177       -4.45 8.25e- 5 1.44e- 4
##  9 ME8        0.0778   -0.0894  0.167        4.24 1.25e- 4 1.95e- 4
## 10 ME3       -0.0774    0.0890 -0.166       -4.04 3.16e- 4 4.42e- 4
## 11 ME11      -0.0737    0.0848 -0.159       -3.77 6.82e- 4 8.68e- 4
## 12 ME4       -0.0606    0.0697 -0.130       -2.92 6.49e- 3 7.19e- 3
## 13 ME2       -0.0597    0.0686 -0.128       -2.89 6.68e- 3 7.19e- 3
## 14 ME12       0.0563   -0.0648  0.121        2.71 1.04e- 2 1.04e- 2

Visualization

Module Eigengene Heatmap

ME_heatmap_data <- ME_data_all %>%
  select(Sample, Genotype, leaf_tissue, Treatment, all_of(module_names_all))

anno_df <- ME_heatmap_data %>%
  select(Sample, Genotype, leaf_tissue, Treatment) %>%
  column_to_rownames("Sample")

ME_matrix <- ME_heatmap_data %>%
  select(all_of(module_names_all)) %>%
  as.matrix()

rownames(ME_matrix) <- ME_heatmap_data$Sample

pheatmap(
  t(ME_matrix),
  annotation_col = anno_df,
  scale = "row",
  clustering_distance_rows = "correlation",
  clustering_distance_cols = "euclidean",
  main = "Module Eigengenes by Sample",
  fontsize = 10
)

Module-Genotype Association

sig_module_plots <- sig_modules %>%
  pull(module) %>%
  head(4)

if (length(sig_module_plots) > 0) {
  plot_data <- ME_data_all %>%
    select(Sample, Genotype, leaf_tissue, all_of(sig_module_plots)) %>%
    pivot_longer(
      cols = all_of(sig_module_plots),
      names_to = "module",
      values_to = "eigengene"
    )
  
  ggplot(plot_data, aes(x = Genotype, y = eigengene, fill = Genotype)) +
    geom_boxplot() +
    geom_jitter(width = 0.2, alpha = 0.5) +
    facet_wrap(~ module, scales = "free_y") +
    theme_classic(base_size = 12) +
    labs(
      title = "Top Genotype-Associated Modules",
      y = "Module Eigengene",
      x = NULL
    ) +
    theme(legend.position = "none")
}

Module Annotation

Combined Network Modules

mod_all <- data.frame(
  gene = names(net_all$colors),
  module = paste0("ME", net_all$colors)
)

module_deg_stats_all <- mod_all %>%
  left_join(
    effects_df %>%
      filter(predictor == "Inv4m") %>%
      select(gene, logFC, adj.P.Val, regulation, in_Inv4m, in_cis),
    by = "gene"
  )

cat("Module assignments (all samples) with DEG statistics:\n")

## Module assignments (all samples) with DEG statistics:

cat("  Total genes:", nrow(module_deg_stats_all), "\n")

##   Total genes: 776

cat("  Genes with stats:", sum(!is.na(module_deg_stats_all$logFC)), "\n")

##   Genes with stats: 776

Genotype-Specific Modules

mod_ctrl <- data.frame(
  gene = names(net_ctrl$colors),
  module = paste0("ME", net_ctrl$colors)
)

module_deg_stats_ctrl <- mod_ctrl %>%
  left_join(
    effects_df %>%
      filter(predictor == "Inv4m") %>%
      select(gene, logFC, adj.P.Val, regulation, in_Inv4m, in_cis),
    by = "gene"
  )

mod_inv4m <- data.frame(
  gene = names(net_inv4m$colors),
  module = paste0("ME", net_inv4m$colors)
)

module_deg_stats_inv4m <- mod_inv4m %>%
  left_join(
    effects_df %>%
      filter(predictor == "Inv4m") %>%
      select(gene, logFC, adj.P.Val, regulation, in_Inv4m, in_cis),
    by = "gene"
  )

cat("\nModule assignments (CTRL) with DEG statistics:\n")

## 
## Module assignments (CTRL) with DEG statistics:

cat("  Total genes:", nrow(module_deg_stats_ctrl), "\n")

##   Total genes: 776

cat("  Genes with stats:", sum(!is.na(module_deg_stats_ctrl$logFC)), "\n")

##   Genes with stats: 776

cat("\nModule assignments (Inv4m) with DEG statistics:\n")

## 
## Module assignments (Inv4m) with DEG statistics:

cat("  Total genes:", nrow(module_deg_stats_inv4m), "\n")

##   Total genes: 776

cat("  Genes with stats:", sum(!is.na(module_deg_stats_inv4m$logFC)), "\n")

##   Genes with stats: 776

Build Edge Lists

cat("\nCalculating adjacency matrix (all samples)...\n")

## 
## Calculating adjacency matrix (all samples)...

adj_all <- adjacency(
  datExpr,
  power = pow_all,
  type = "unsigned"
)

cat("Adjacency matrix dimensions:", dim(adj_all), "\n")

## Adjacency matrix dimensions: 776 776

edge_indices_all <- which(adj_all > 0 & upper.tri(adj_all, diag = FALSE), 
                          arr.ind = TRUE)

edge_list_all <- data.frame(
  from = rownames(adj_all)[edge_indices_all[, 1]],
  to = colnames(adj_all)[edge_indices_all[, 2]],
  weight = adj_all[edge_indices_all]
) %>%
  arrange(desc(weight))

cat("Total edges (weight > 0):", nrow(edge_list_all), "\n")

## Total edges (weight > 0): 300700

cat("\nCalculating adjacency matrices (genotype-specific)...\n")

## 
## Calculating adjacency matrices (genotype-specific)...

adj_ctrl <- adjacency(
  datExpr_ctrl,
  power = pow_ctrl,
  type = "unsigned"
)

adj_inv4m <- adjacency(
  datExpr_inv4m,
  power = pow_inv4m,
  type = "unsigned"
)

cat("CTRL adjacency dimensions:", dim(adj_ctrl), "\n")

## CTRL adjacency dimensions: 776 776

cat("Inv4m adjacency dimensions:", dim(adj_inv4m), "\n")

## Inv4m adjacency dimensions: 776 776

Module Summary Statistics

Combined Network

modules_all_summary <- module_deg_stats_all %>%
  group_by(module) %>%
  summarise(
    n_genes = n(),
    n_upregulated = sum(regulation == "Upregulated", na.rm = TRUE),
    n_downregulated = sum(regulation == "Downregulated", na.rm = TRUE),
    mean_logFC = mean(logFC, na.rm = TRUE),
    median_logFC = median(logFC, na.rm = TRUE),
    n_in_inv4m = sum(in_Inv4m, na.rm = TRUE),
    n_in_cis = sum(in_cis, na.rm = TRUE),
    pct_in_cis = 100 * n_in_cis / n_genes,
    .groups = "drop"
  ) %>%
  arrange(desc(abs(mean_logFC)))

cat("Module summary (all samples):\n")

## Module summary (all samples):

print(modules_all_summary)

## # A tibble: 14 × 9
##    module n_genes n_upregulated n_downregulated mean_logFC median_logFC
##    <chr>    <int>         <int>           <int>      <dbl>        <dbl>
##  1 ME14        10             1               9   -0.711         -0.370
##  2 ME8         30            28               2    0.588          0.548
##  3 ME9         29             8              21   -0.430         -0.434
##  4 ME7         31             0              31   -0.388         -0.325
##  5 ME10        29             0              29   -0.363         -0.294
##  6 ME3         60             4              56   -0.319         -0.287
##  7 ME2         77             5              72   -0.300         -0.373
##  8 ME4         50             4              46   -0.221         -0.316
##  9 ME1        340           140             200   -0.189         -0.350
## 10 ME12        20            11               9    0.182          0.383
## 11 ME11        21             5              16   -0.168         -0.321
## 12 ME13        13             6               7   -0.155         -0.279
## 13 ME5         35             9              26   -0.155         -0.410
## 14 ME6         31            16              15   -0.00166        0.335
## # ℹ 3 more variables: n_in_inv4m <int>, n_in_cis <int>, pct_in_cis <dbl>

write_csv(modules_all_summary, file.path(paths$intermediate, "modules_all_summary.csv"))

Genotype-Specific Networks

modules_ctrl_summary <- module_deg_stats_ctrl %>%
  group_by(module) %>%
  summarise(
    n_genes = n(),
    n_upregulated = sum(regulation == "Upregulated", na.rm = TRUE),
    n_downregulated = sum(regulation == "Downregulated", na.rm = TRUE),
    mean_logFC = mean(logFC, na.rm = TRUE),
    median_logFC = median(logFC, na.rm = TRUE),
    n_in_inv4m = sum(in_Inv4m, na.rm = TRUE),
    n_in_cis = sum(in_cis, na.rm = TRUE),
    pct_in_cis = 100 * n_in_cis / n_genes,
    .groups = "drop"
  ) %>%
  arrange(desc(abs(mean_logFC)))

modules_inv4m_summary <- module_deg_stats_inv4m %>%
  group_by(module) %>%
  summarise(
    n_genes = n(),
    n_upregulated = sum(regulation == "Upregulated", na.rm = TRUE),
    n_downregulated = sum(regulation == "Downregulated", na.rm = TRUE),
    mean_logFC = mean(logFC, na.rm = TRUE),
    median_logFC = median(logFC, na.rm = TRUE),
    n_in_inv4m = sum(in_Inv4m, na.rm = TRUE),
    n_in_cis = sum(in_cis, na.rm = TRUE),
    pct_in_cis = 100 * n_in_cis / n_genes,
    .groups = "drop"
  ) %>%
  arrange(desc(abs(mean_logFC)))

cat("\nModule summary (CTRL):\n")

## 
## Module summary (CTRL):

print(modules_ctrl_summary)

## # A tibble: 19 × 9
##    module n_genes n_upregulated n_downregulated mean_logFC median_logFC
##    <chr>    <int>         <int>           <int>      <dbl>        <dbl>
##  1 ME16        15            12               3     3.21          4.98 
##  2 ME0         27            17              10     0.984         0.763
##  3 ME11        28             7              21    -0.857        -0.399
##  4 ME8         36            12              24    -0.819        -0.428
##  5 ME17        15             4              11     0.752        -0.589
##  6 ME9         36            16              20    -0.702        -0.360
##  7 ME10        30             4              26    -0.683        -0.326
##  8 ME2        102            27              75    -0.655        -0.351
##  9 ME7         39            10              29    -0.524        -0.265
## 10 ME1        116            13             103    -0.509        -0.365
## 11 ME5         56            27              29     0.426        -0.231
## 12 ME18        13             2              11    -0.245        -0.305
## 13 ME4         63            36              27     0.236         0.478
## 14 ME6         40            18              22    -0.178        -0.501
## 15 ME3         82            13              69    -0.173        -0.305
## 16 ME12        27             5              22    -0.115        -0.242
## 17 ME15        16             3              13    -0.0859       -0.376
## 18 ME13        18             7              11     0.0708       -0.251
## 19 ME14        17             4              13    -0.0506       -0.454
## # ℹ 3 more variables: n_in_inv4m <int>, n_in_cis <int>, pct_in_cis <dbl>

cat("\nModule summary (Inv4m):\n")

## 
## Module summary (Inv4m):

print(modules_inv4m_summary)

## # A tibble: 19 × 9
##    module n_genes n_upregulated n_downregulated mean_logFC median_logFC
##    <chr>    <int>         <int>           <int>      <dbl>        <dbl>
##  1 ME12        16             1              15    -6.32         -6.73 
##  2 ME8         26            22               4     1.93          0.730
##  3 ME7         33            18              15     1.50          0.272
##  4 ME18        11            10               1     0.923         0.858
##  5 ME11        17            10               7     0.862         0.534
##  6 ME13        15             0              15    -0.839        -0.312
##  7 ME17        12            10               2     0.689         0.686
##  8 ME6         38             1              37    -0.687        -0.390
##  9 ME16        12             3               9    -0.416        -0.662
## 10 ME10        19             1              18    -0.416        -0.323
## 11 ME3         45             1              44    -0.404        -0.294
## 12 ME0        288            89             199    -0.352        -0.370
## 13 ME2         50            22              28     0.300        -0.199
## 14 ME9         21             3              18    -0.241        -0.374
## 15 ME5         40            13              27    -0.209        -0.361
## 16 ME1         62            14              48    -0.201        -0.315
## 17 ME15        13             4               9     0.155        -0.263
## 18 ME14        14             2              12    -0.0946       -0.350
## 19 ME4         44            13              31     0.0120       -0.390
## # ℹ 3 more variables: n_in_inv4m <int>, n_in_cis <int>, pct_in_cis <dbl>

write_csv(modules_ctrl_summary, file.path(paths$intermediate, "modules_ctrl_summary.csv"))
write_csv(modules_inv4m_summary, file.path(paths$intermediate, "modules_inv4m_summary.csv"))

Identify Hub Genes

Combined Network

cat("\nCalculating intramodular connectivity (all samples)...\n")

## 
## Calculating intramodular connectivity (all samples)...

connectivity_all <- intramodularConnectivity(adj_all, net_all$colors)

hubs_all <- mod_all %>%
  mutate(
    module_numeric = as.numeric(gsub("ME", "", module))
  ) %>%
  left_join(
    connectivity_all %>%
      rownames_to_column("gene") %>%
      as_tibble(),
    by = "gene"
  ) %>%
  left_join(
    effects_df %>%
      filter(predictor == "Inv4m") %>%
      select(gene, locus_label, desc_merged, logFC, adj.P.Val, 
             regulation, in_Inv4m, in_cis),
    by = "gene"
  ) %>%
  arrange(module, desc(kWithin))

cat("Hub genes calculated (all samples)\n")

## Hub genes calculated (all samples)

cat("  Mean kWithin:", round(mean(hubs_all$kWithin, na.rm = TRUE), 3), "\n")

##   Mean kWithin: 15.734

cat("  Median kWithin:", round(median(hubs_all$kWithin, na.rm = TRUE), 3), "\n")

##   Median kWithin: 6.15

hubs_all_top <- hubs_all %>%
  filter(abs(logFC) >= 1.5) %>%
  group_by(module) %>%
  arrange(desc(kWithin)) %>%
  slice_head(n = 10) %>%
  ungroup()

cat("\nTop 10 hub genes per module (all samples):\n")

## 
## Top 10 hub genes per module (all samples):

print(hubs_all_top %>% 
        select(module, gene, locus_label, kWithin, logFC, in_cis))

## # A tibble: 21 × 6
##    module gene            locus_label kWithin logFC in_cis
##    <chr>  <chr>           <chr>         <dbl> <dbl> <lgl> 
##  1 ME1    Zm00001eb264460 ychf           89.8  6.54 FALSE 
##  2 ME1    Zm00001eb213070 gbp1           89.5  8.97 FALSE 
##  3 ME1    Zm00001eb189640 fam32a_2       87.6  5.58 TRUE  
##  4 ME1    Zm00001eb192880 tet2           86.8 -6.95 TRUE  
##  5 ME1    Zm00001eb192850 <NA>           86.4 -1.99 TRUE  
##  6 ME1    Zm00001eb187280 gpm458         85.7  9.91 TRUE  
##  7 ME1    Zm00001eb196760 <NA>           85.0 -8.51 TRUE  
##  8 ME1    Zm00001eb126050 pip5k1         84.9  4.71 FALSE 
##  9 ME1    Zm00001eb190240 <NA>           83.6 -9.33 TRUE  
## 10 ME1    Zm00001eb191790 jmj2           82.4 -3.55 TRUE  
## # ℹ 11 more rows

write_csv(hubs_all, file.path(paths$intermediate, "hubs_all.csv"))
write_csv(hubs_all_top, file.path(paths$intermediate, "hubs_all_top.csv"))

cat("\nHub gene files exported (all samples)\n")

## 
## Hub gene files exported (all samples)

connectivity_threshold_all <- quantile(hubs_all$kWithin, 0.75, na.rm = TRUE)

hubs_all_cis <- hubs_all %>%
  filter(in_cis == TRUE, kWithin > connectivity_threshold_all) %>%
  arrange(desc(kWithin))

cat("\nCis hub genes (all samples, kWithin > 75th percentile =", 
    round(connectivity_threshold_all, 2), "):", nrow(hubs_all_cis), "\n")

## 
## Cis hub genes (all samples, kWithin > 75th percentile = 16.26 ): 149

if (nrow(hubs_all_cis) > 0) {
  print(hubs_all_cis %>% 
          select(module, gene, locus_label, kWithin, logFC, in_Inv4m))
  
  write_csv(hubs_all_cis, file.path(paths$intermediate, "hubs_all_cis.csv"))
}

##     module            gene locus_label  kWithin      logFC in_Inv4m
## 1      ME1 Zm00001eb189640    fam32a_2 87.60824  5.5840769    FALSE
## 2      ME1 Zm00001eb192880        tet2 86.82678 -6.9457261     TRUE
## 3      ME1 Zm00001eb192850        <NA> 86.39492 -1.9946674     TRUE
## 4      ME1 Zm00001eb187280      gpm458 85.67797  9.9097830    FALSE
## 5      ME1 Zm00001eb196760        <NA> 85.01359 -8.5057276    FALSE
## 6      ME1 Zm00001eb190240        <NA> 83.58921 -9.3262051    FALSE
## 7      ME1 Zm00001eb191790        jmj2 82.43619 -3.5479196     TRUE
## 8      ME1 Zm00001eb191820        jmj4 81.50158 -2.8826897     TRUE
## 9      ME1 Zm00001eb192630        <NA> 80.50286 -1.7229795     TRUE
## 10     ME1 Zm00001eb188750    hsp70-17 79.53606 -6.5893888    FALSE
## 11     ME1 Zm00001eb189580        <NA> 78.83744  7.5068431    FALSE
## 12     ME1 Zm00001eb189350        <NA> 78.79253 -9.6708956    FALSE
## 13     ME1 Zm00001eb188420        <NA> 78.02280 -6.8786411    FALSE
## 14     ME1 Zm00001eb188430        <NA> 78.02280 -6.8786411    FALSE
## 15     ME1 Zm00001eb190750       jmj21 77.44843  2.3124566     TRUE
## 16     ME1 Zm00001eb193870      ftshi3 77.25884 -5.4979858     TRUE
## 17     ME1 Zm00001eb193260       trxl2 77.15654 -6.0252641     TRUE
## 18     ME1 Zm00001eb188240    ropgef14 75.78290 -3.0530829    FALSE
## 19     ME1 Zm00001eb190670       zfrvt 75.23342  1.9515558     TRUE
## 20     ME1 Zm00001eb192170        <NA> 75.18992 -7.4688231     TRUE
## 21     ME1 Zm00001eb192970        <NA> 74.48943 -6.3077016     TRUE
## 22     ME1 Zm00001eb192330      roc4-1 73.64992  5.0518680     TRUE
## 23     ME1 Zm00001eb192080        <NA> 73.55259 -1.4839851     TRUE
## 24     ME1 Zm00001eb188070        arf1 73.20311 -2.8277808    FALSE
## 25     ME1 Zm00001eb188460       dof43 72.69711 -5.1335736    FALSE
## 26     ME1 Zm00001eb189200        <NA> 72.62853 -6.1597046    FALSE
## 27     ME1 Zm00001eb189170        <NA> 72.07607 -1.8651869    FALSE
## 28     ME1 Zm00001eb192160        <NA> 72.06273 -6.8776969     TRUE
## 29     ME1 Zm00001eb190350        <NA> 71.92532 -1.8752894    FALSE
## 30     ME1 Zm00001eb193620      zfcchc 71.71367 -8.8405439     TRUE
## 31     ME1 Zm00001eb193960        <NA> 70.54737 -1.2610279     TRUE
## 32     ME1 Zm00001eb193130        <NA> 70.09817 -5.6331437     TRUE
## 33     ME1 Zm00001eb188610        <NA> 69.58071 -5.9771086    FALSE
## 34     ME1 Zm00001eb194200        <NA> 69.01066 -1.3199224     TRUE
## 35     ME1 Zm00001eb188570      his2b5 67.77224 -5.3380982    FALSE
## 36     ME1 Zm00001eb189910      metrs1 67.49021 -2.9121795    FALSE
## 37     ME1 Zm00001eb195370        <NA> 67.29793  1.6383226    FALSE
## 38     ME1 Zm00001eb188330        <NA> 66.21055 -1.1528583    FALSE
## 39     ME1 Zm00001eb189820        <NA> 66.13757  5.2973522    FALSE
## 40     ME1 Zm00001eb195120        <NA> 63.96160 -1.2011818    FALSE
## 41     ME1 Zm00001eb190090        <NA> 63.91537  7.1778814    FALSE
## 42     ME1 Zm00001eb189570        <NA> 63.76418  5.5151847    FALSE
## 43     ME1 Zm00001eb192180        <NA> 63.52667 -5.1958191     TRUE
## 44     ME1 Zm00001eb196780        <NA> 63.50148 -6.8610167    FALSE
## 45     ME1 Zm00001eb195860        <NA> 62.97926 -1.4718951    FALSE
## 46     ME1 Zm00001eb192980        <NA> 62.90425 -1.1956236     TRUE
## 47     ME1 Zm00001eb187430      his402 62.49404 -8.1694502    FALSE
## 48     ME1 Zm00001eb196600        r3h1 62.36288 -2.2746235    FALSE
## 49     ME1 Zm00001eb195800        <NA> 61.61740  1.6858934    FALSE
## 50     ME1 Zm00001eb196830        <NA> 58.89420 -5.3301933    FALSE
## 51     ME1 Zm00001eb188220        <NA> 57.93621 -2.9845203    FALSE
## 52     ME1 Zm00001eb191990        pig3 56.96594  1.9529439     TRUE
## 53     ME1 Zm00001eb192900        <NA> 56.82845 -0.8713048     TRUE
## 54     ME1 Zm00001eb187890         smn 56.73194  3.9673651    FALSE
## 55     ME1 Zm00001eb189900       hb115 56.64169 -2.0719433    FALSE
## 56     ME1 Zm00001eb188990        fad1 55.83089 -2.1286584    FALSE
## 57     ME1 Zm00001eb191260        <NA> 55.65987 -1.3868319     TRUE
## 58     ME1 Zm00001eb191320        <NA> 55.17959 -1.2326403     TRUE
## 59     ME1 Zm00001eb194300       paspa 54.89177 -4.0237443     TRUE
## 60     ME1 Zm00001eb194700        <NA> 54.22765 -2.1314797     TRUE
## 61     ME1 Zm00001eb195080        <NA> 54.10796 -4.8397541    FALSE
## 62     ME1 Zm00001eb196200        <NA> 53.91908 -2.3537716    FALSE
## 63     ME1 Zm00001eb196170       mfsd2 53.81298 -1.6469160    FALSE
## 64     ME1 Zm00001eb194380        sec6 53.67194  2.8225461     TRUE
## 65     ME1 Zm00001eb195940        <NA> 53.59613 -0.6912339    FALSE
## 66     ME1 Zm00001eb189390        <NA> 52.86676 -1.4479467    FALSE
## 67     ME1 Zm00001eb195560        <NA> 52.28661 -5.7425424    FALSE
## 68     ME1 Zm00001eb191710        <NA> 52.07647 -0.7872265     TRUE
## 69     ME1 Zm00001eb196880      npi270 51.14548  1.3003241    FALSE
## 70     ME1 Zm00001eb193300        <NA> 50.09204 -4.8842457     TRUE
## 71     ME1 Zm00001eb189700        <NA> 49.47098  4.4551791    FALSE
## 72     ME1 Zm00001eb190500        <NA> 48.73563 -0.5672530     TRUE
## 73     ME1 Zm00001eb189420        <NA> 48.51522 -0.6316545    FALSE
## 74     ME1 Zm00001eb196100        abf2 48.35790  1.7886406    FALSE
## 75     ME1 Zm00001eb196920        <NA> 48.15056 -2.8468789    FALSE
## 76     ME1 Zm00001eb187030        <NA> 47.89935 -3.2839911    FALSE
## 77     ME1 Zm00001eb195790        <NA> 46.91848  3.1218159    FALSE
## 78     ME1 Zm00001eb196240        pin4 46.62636 -3.3660865    FALSE
## 79     ME1 Zm00001eb191430        <NA> 45.97349 -1.0298254     TRUE
## 80     ME1 Zm00001eb189920       aldh2 45.25686  2.6863326    FALSE
## 81     ME1 Zm00001eb192120        <NA> 44.43765 -2.3947086     TRUE
## 82     ME1 Zm00001eb191150        <NA> 44.17096  1.4618638     TRUE
## 83     ME1 Zm00001eb192250        <NA> 44.16916 -0.8130420     TRUE
## 84     ME1 Zm00001eb188030        <NA> 43.46637  4.2180694    FALSE
## 85     ME1 Zm00001eb195970        <NA> 43.25589 -2.6358333    FALSE
## 86     ME1 Zm00001eb189710        <NA> 43.25525  3.7158176    FALSE
## 87     ME1 Zm00001eb193120        <NA> 42.19165  2.3770809     TRUE
## 88     ME1 Zm00001eb187870        <NA> 40.95528  2.5341132    FALSE
## 89     ME1 Zm00001eb195900        ms30 40.14422 -2.8851397    FALSE
## 90     ME1 Zm00001eb196090      pmei47 40.08216  5.3555573    FALSE
## 91     ME1 Zm00001eb191620        <NA> 39.12616 -1.4519423     TRUE
## 92     ME1 Zm00001eb189120        <NA> 39.01994 -0.5983260    FALSE
## 93     ME1 Zm00001eb192580        <NA> 38.89437  2.7379520     TRUE
## 94     ME1 Zm00001eb195600        ptk2 38.69780 -1.6043996    FALSE
## 95     ME1 Zm00001eb190660        <NA> 38.13736  0.8335129     TRUE
## 96     ME1 Zm00001eb192360        <NA> 37.18787 -0.8833142     TRUE
## 97     ME1 Zm00001eb194180        <NA> 35.83976  4.0975260     TRUE
## 98     ME1 Zm00001eb188700     hsp70-7 35.55379 -1.0096395    FALSE
## 99     ME1 Zm00001eb193470        <NA> 34.95628  1.4143723     TRUE
## 100    ME1 Zm00001eb191560        <NA> 34.69698 -1.3020994     TRUE
## 101    ME1 Zm00001eb197020        <NA> 34.14726 -0.9050130    FALSE
## 102    ME1 Zm00001eb192350        <NA> 33.58008 -2.2712141     TRUE
## 103    ME1 Zm00001eb188680        <NA> 33.08431 -0.8451098    FALSE
## 104    ME1 Zm00001eb190880        <NA> 32.27555  0.4458239     TRUE
## 105    ME1 Zm00001eb193560        <NA> 32.26116 -0.4043503     TRUE
## 106    ME1 Zm00001eb186400        <NA> 31.94061  7.7457413    FALSE
## 107    ME1 Zm00001eb196500        <NA> 31.75273 -1.4366628    FALSE
## 108    ME1 Zm00001eb192460        <NA> 31.57491  0.8579075     TRUE
## 109    ME1 Zm00001eb192400        <NA> 30.74269 -0.8670134     TRUE
## 110    ME1 Zm00001eb195050        <NA> 30.11611 -0.3595747    FALSE
## 111    ME1 Zm00001eb195620        <NA> 29.69879 -1.8772725    FALSE
## 112    ME1 Zm00001eb196180        tom2 29.66175 -1.6691084    FALSE
## 113    ME1 Zm00001eb193970       fbl41 29.47838 -1.2902767     TRUE
## 114    ME1 Zm00001eb195290        <NA> 28.37934  1.4354492    FALSE
## 115    ME1 Zm00001eb186390        <NA> 27.82918  5.2882507    FALSE
## 116    ME1 Zm00001eb194240        <NA> 26.86493  0.7327490     TRUE
## 117    ME1 Zm00001eb195830       zhd18 26.85401 -4.3885577    FALSE
## 118    ME1 Zm00001eb186410        <NA> 26.23715  4.9770347    FALSE
## 119    ME1 Zm00001eb196550        <NA> 26.14473 -1.3159081    FALSE
## 120    ME1 Zm00001eb194270       flz22 26.00829  2.6619638     TRUE
## 121    ME1 Zm00001eb190190        <NA> 25.87007  0.5405926    FALSE
## 122    ME1 Zm00001eb195480        <NA> 25.73144 -0.2307286    FALSE
## 123    ME1 Zm00001eb192840        <NA> 25.10500 -1.6572585     TRUE
## 124    ME1 Zm00001eb196630       ufgt4 24.83341  3.9050799    FALSE
## 125    ME1 Zm00001eb189510        dof9 24.39562  0.5734774    FALSE
## 126    ME1 Zm00001eb195180        <NA> 24.29245 -2.3220037    FALSE
## 127    ME1 Zm00001eb188890        <NA> 23.53009 -0.5838285    FALSE
## 128    ME1 Zm00001eb192060        <NA> 23.43976 -0.3905260     TRUE
## 129    ME1 Zm00001eb197300        <NA> 23.13444 -2.6729304    FALSE
## 130    ME1 Zm00001eb193510        <NA> 21.54209 -0.5357527     TRUE
## 131    ME1 Zm00001eb194590        <NA> 20.65390 -0.8451295     TRUE
## 132    ME1 Zm00001eb191680        <NA> 20.39198 -0.3345883     TRUE
## 133    ME1 Zm00001eb188550 rz567b(klc) 20.23608  2.3169023    FALSE
## 134    ME1 Zm00001eb194080        <NA> 20.08124 -4.9094407     TRUE
## 135    ME1 Zm00001eb194610        <NA> 19.37911 -2.7287497     TRUE
## 136    ME1 Zm00001eb195960        <NA> 19.36948  0.7766816    FALSE
## 137    ME1 Zm00001eb192340        vps3 18.61204 -0.6099370     TRUE
## 138    ME1 Zm00001eb194570        <NA> 18.57230 -1.7845376     TRUE
## 139    ME1 Zm00001eb193660        nii2 18.30599 -1.9515935     TRUE
## 140    ME1 Zm00001eb186450        <NA> 18.15185 -0.3732141    FALSE
## 141    ME1 Zm00001eb187460        <NA> 17.96826  2.1857894    FALSE
## 142    ME1 Zm00001eb191780     nactf51 17.67973 -3.5605202     TRUE
## 143    ME1 Zm00001eb188050        <NA> 17.42347 -1.3477292    FALSE
## 144    ME1 Zm00001eb190450        <NA> 17.06172 -1.5667953    FALSE
## 145    ME1 Zm00001eb194600        <NA> 17.02625 -0.8923818     TRUE
## 146    ME1 Zm00001eb191000        <NA> 16.80309 -3.5079674     TRUE
## 147    ME1 Zm00001eb188350        <NA> 16.72521 -0.8678902    FALSE
## 148    ME1 Zm00001eb192260        <NA> 16.55576 -0.4354540     TRUE
## 149    ME1 Zm00001eb194680        rop8 16.38496 -1.6676644     TRUE

Genotype-Specific Networks

cat("\nCalculating intramodular connectivity (genotype-specific)...\n")

## 
## Calculating intramodular connectivity (genotype-specific)...

connectivity_ctrl <- intramodularConnectivity(adj_ctrl, net_ctrl$colors)
connectivity_inv4m <- intramodularConnectivity(adj_inv4m, net_inv4m$colors)

hubs_ctrl <- mod_ctrl %>%
  mutate(
    module_numeric = as.numeric(gsub("ME", "", module))
  ) %>%
  left_join(
    connectivity_ctrl %>%
      rownames_to_column("gene") %>%
      as_tibble(),
    by = "gene"
  ) %>%
  left_join(
    effects_df %>%
      filter(predictor == "Inv4m") %>%
      select(gene, locus_label, desc_merged, logFC, adj.P.Val, 
             regulation, in_Inv4m, in_cis),
    by = "gene"
  ) %>%
  arrange(module, desc(kWithin))

hubs_inv4m <- mod_inv4m %>%
  mutate(
    module_numeric = as.numeric(gsub("ME", "", module))
  ) %>%
  left_join(
    connectivity_inv4m %>%
      rownames_to_column("gene") %>%
      as_tibble(),
    by = "gene"
  ) %>%
  left_join(
    effects_df %>%
      filter(predictor == "Inv4m") %>%
      select(gene, locus_label, desc_merged, logFC, adj.P.Val, 
             regulation, in_Inv4m, in_cis),
    by = "gene"
  ) %>%
  arrange(module, desc(kWithin))

cat("Hub genes calculated\n")

## Hub genes calculated

cat("  CTRL - Mean kWithin:", round(mean(hubs_ctrl$kWithin, na.rm = TRUE), 3), "\n")

##   CTRL - Mean kWithin: 6.166

cat("  CTRL - Median kWithin:", round(median(hubs_ctrl$kWithin, na.rm = TRUE), 3), "\n")

##   CTRL - Median kWithin: 3.742

cat("  Inv4m - Mean kWithin:", round(mean(hubs_inv4m$kWithin, na.rm = TRUE), 3), "\n")

##   Inv4m - Mean kWithin: 0.443

cat("  Inv4m - Median kWithin:", round(median(hubs_inv4m$kWithin, na.rm = TRUE), 3), "\n")

##   Inv4m - Median kWithin: 0.164

hubs_ctrl_top <- hubs_ctrl %>%
  filter(abs(logFC) >= 1.5) %>%
  group_by(module) %>%
  arrange(desc(kWithin)) %>%
  slice_head(n = 10) %>%
  ungroup()

hubs_inv4m_top <- hubs_inv4m %>%
  filter(abs(logFC) >= 1.5) %>%
  group_by(module) %>%
  arrange(desc(kWithin)) %>%
  slice_head(n = 10) %>%
  ungroup()

cat("\nTop 10 hub genes per module (CTRL):\n")

## 
## Top 10 hub genes per module (CTRL):

print(hubs_ctrl_top %>% 
        select(module, gene, locus_label, kWithin, logFC, in_cis))

## # A tibble: 123 × 6
##    module gene            locus_label kWithin logFC in_cis
##    <chr>  <chr>           <chr>         <dbl> <dbl> <lgl> 
##  1 ME0    Zm00001eb187440 <NA>         0.174  -2.00 TRUE  
##  2 ME0    Zm00001eb166340 <NA>         0.154  -4.62 FALSE 
##  3 ME0    Zm00001eb192580 <NA>         0.113   2.74 TRUE  
##  4 ME0    Zm00001eb189820 <NA>         0.112   5.30 TRUE  
##  5 ME0    Zm00001eb187870 <NA>         0.107   2.53 TRUE  
##  6 ME0    Zm00001eb334460 <NA>         0.0673  4.25 FALSE 
##  7 ME0    Zm00001eb186840 <NA>         0.0592  1.69 TRUE  
##  8 ME0    Zm00001eb189700 <NA>         0.0544  4.46 TRUE  
##  9 ME0    Zm00001eb194570 <NA>         0.0419 -1.78 TRUE  
## 10 ME0    Zm00001eb257900 mcm5         0.0395 -2.74 FALSE 
## # ℹ 113 more rows

cat("\nTop 10 hub genes per module (Inv4m):\n")

## 
## Top 10 hub genes per module (Inv4m):

print(hubs_inv4m_top %>% 
        select(module, gene, locus_label, kWithin, logFC, in_cis))

## # A tibble: 88 × 6
##    module gene            locus_label kWithin logFC in_cis
##    <chr>  <chr>           <chr>         <dbl> <dbl> <lgl> 
##  1 ME0    Zm00001eb334460 <NA>         0.0821  4.25 FALSE 
##  2 ME0    Zm00001eb371380 <NA>         0.0613  3.36 FALSE 
##  3 ME0    Zm00001eb157030 <NA>         0.0577  5.74 FALSE 
##  4 ME0    Zm00001eb195370 <NA>         0.0438  1.64 TRUE  
##  5 ME0    Zm00001eb197300 <NA>         0.0423 -2.67 TRUE  
##  6 ME0    Zm00001eb195620 <NA>         0.0373 -1.88 TRUE  
##  7 ME0    Zm00001eb041890 kan1         0.0365  1.61 FALSE 
##  8 ME0    Zm00001eb190750 jmj21        0.0296  2.31 TRUE  
##  9 ME0    Zm00001eb190450 <NA>         0.0266 -1.57 TRUE  
## 10 ME0    Zm00001eb248260 <NA>         0.0264  2.08 FALSE 
## # ℹ 78 more rows

write_csv(hubs_ctrl, file.path(paths$intermediate, "hubs_ctrl.csv"))
write_csv(hubs_ctrl_top, file.path(paths$intermediate, "hubs_ctrl_top.csv"))
write_csv(hubs_inv4m, file.path(paths$intermediate, "hubs_inv4m.csv"))
write_csv(hubs_inv4m_top, file.path(paths$intermediate, "hubs_inv4m_top.csv"))

cat("\nHub gene files exported (genotype-specific)\n")

## 
## Hub gene files exported (genotype-specific)

connectivity_threshold_ctrl <- quantile(hubs_ctrl$kWithin, 0.75, na.rm = TRUE)
connectivity_threshold_inv4m <- quantile(hubs_inv4m$kWithin, 0.75, na.rm = TRUE)

hubs_ctrl_cis <- hubs_ctrl %>%
  filter(in_cis == TRUE, kWithin > connectivity_threshold_ctrl) %>%
  arrange(desc(kWithin))

hubs_inv4m_cis <- hubs_inv4m %>%
  filter(in_cis == TRUE, kWithin > connectivity_threshold_inv4m) %>%
  arrange(desc(kWithin))

cat("\nCis hub genes (CTRL, kWithin > 75th percentile =", 
    round(connectivity_threshold_ctrl, 2), "):", nrow(hubs_ctrl_cis), "\n")

## 
## Cis hub genes (CTRL, kWithin > 75th percentile = 7.68 ): 52

if (nrow(hubs_ctrl_cis) > 0) {
  print(hubs_ctrl_cis %>% 
          select(module, gene, locus_label, kWithin, logFC, in_Inv4m))
  
  write_csv(hubs_ctrl_cis, file.path(paths$intermediate, "hubs_ctrl_cis.csv"))
}

##    module            gene locus_label   kWithin      logFC in_Inv4m
## 1     ME1 Zm00001eb192310        <NA> 25.475435 -0.4493898     TRUE
## 2     ME1 Zm00001eb192140        <NA> 22.674525 -0.5957101     TRUE
## 3     ME2 Zm00001eb191790        jmj2 19.354256 -3.5479196     TRUE
## 4     ME3 Zm00001eb195970        <NA> 16.699627 -2.6358333    FALSE
## 5     ME2 Zm00001eb195960        <NA> 15.448279  0.7766816    FALSE
## 6     ME2 Zm00001eb196080        <NA> 15.390397  0.4392978    FALSE
## 7     ME2 Zm00001eb191820        jmj4 15.355796 -2.8826897     TRUE
## 8     ME3 Zm00001eb189030       tola1 15.270774 -0.4167695    FALSE
## 9     ME5 Zm00001eb190360        <NA> 13.943297 -0.9146507    FALSE
## 10    ME5 Zm00001eb194130       rpl29 13.522970 -0.5654438     TRUE
## 11    ME5 Zm00001eb186900       pip1d 13.258877  1.0867126    FALSE
## 12    ME1 Zm00001eb188520        <NA> 13.028288  0.4027907    FALSE
## 13    ME5 Zm00001eb187200        <NA> 12.951922  0.5562427    FALSE
## 14    ME5 Zm00001eb188310        <NA> 12.566939  0.9907331    FALSE
## 15    ME5 Zm00001eb195780    AY104784 12.443780 -0.4335815    FALSE
## 16    ME5 Zm00001eb196030      cipk33 12.380820 -0.7361484    FALSE
## 17    ME2 Zm00001eb191890         ss5 12.354537 -0.6662528     TRUE
## 18    ME2 Zm00001eb195620        <NA> 12.108495 -1.8772725    FALSE
## 19    ME5 Zm00001eb187590        <NA> 11.750940 -1.4486884    FALSE
## 20    ME5 Zm00001eb189210        <NA> 11.533746 -0.3135679    FALSE
## 21    ME6 Zm00001eb193820        <NA> 11.153724  0.4699150     TRUE
## 22    ME2 Zm00001eb194670        <NA> 10.885101 -0.3568029     TRUE
## 23    ME6 Zm00001eb189950        <NA> 10.834141  0.9848448    FALSE
## 24    ME2 Zm00001eb195120        <NA> 10.558953 -1.2011818    FALSE
## 25    ME5 Zm00001eb188220        <NA> 10.199101 -2.9845203    FALSE
## 26    ME5 Zm00001eb189850        <NA> 10.049028 -0.4890672    FALSE
## 27    ME5 Zm00001eb197120        <NA>  9.980317 -0.8027816    FALSE
## 28    ME2 Zm00001eb192360        <NA>  9.853863 -0.8833142     TRUE
## 29    ME6 Zm00001eb192910     lkrsdh1  9.794229 -1.8177127     TRUE
## 30    ME5 Zm00001eb190390        <NA>  9.705346 -0.2950611    FALSE
## 31    ME2 Zm00001eb196310       emp11  9.580118 -0.5777013    FALSE
## 32    ME8 Zm00001eb189140        <NA>  9.337569 -1.1740956    FALSE
## 33    ME8 Zm00001eb190110        <NA>  9.336529 -0.3443905    FALSE
## 34    ME5 Zm00001eb192110        <NA>  9.186318  0.3908126     TRUE
## 35    ME5 Zm00001eb196390        <NA>  9.157441  0.4858115    FALSE
## 36    ME6 Zm00001eb193660        nii2  9.120091 -1.9515935     TRUE
## 37    ME5 Zm00001eb194270       flz22  8.900508  2.6619638     TRUE
## 38    ME4 Zm00001eb193790         tu1  8.812693 -0.4217700     TRUE
## 39    ME6 Zm00001eb189960        <NA>  8.793860  0.9247815    FALSE
## 40    ME5 Zm00001eb187850        <NA>  8.625734 -3.9035202    FALSE
## 41    ME5 Zm00001eb190190        <NA>  8.614837  0.5405926    FALSE
## 42    ME6 Zm00001eb194560        <NA>  8.391949  0.6795678     TRUE
## 43    ME4 Zm00001eb187910        <NA>  8.302450 -0.9566888    FALSE
## 44    ME6 Zm00001eb196880      npi270  8.208947  1.3003241    FALSE
## 45    ME6 Zm00001eb195600        ptk2  8.104428 -1.6043996    FALSE
## 46    ME8 Zm00001eb193040       cpn10  8.018739  0.6148654     TRUE
## 47    ME8 Zm00001eb192960      ereb14  7.988048 -0.4713370     TRUE
## 48    ME5 Zm00001eb189440        <NA>  7.960522 -2.2523055    FALSE
## 49    ME8 Zm00001eb195380        <NA>  7.956684  0.3763682    FALSE
## 50    ME2 Zm00001eb187970        wtf4  7.870608 -0.6071448    FALSE
## 51    ME6 Zm00001eb192640        <NA>  7.784881  0.7591437     TRUE
## 52    ME5 Zm00001eb187240        <NA>  7.741272 -0.2770947    FALSE

cat("\nCis hub genes (Inv4m, kWithin > 75th percentile =", 
    round(connectivity_threshold_inv4m, 2), "):", nrow(hubs_inv4m_cis), "\n")

## 
## Cis hub genes (Inv4m, kWithin > 75th percentile = 0.57 ): 75

if (nrow(hubs_inv4m_cis) > 0) {
  print(hubs_inv4m_cis %>% 
          select(module, gene, locus_label, kWithin, logFC, in_Inv4m))
  
  write_csv(hubs_inv4m_cis, file.path(paths$intermediate, "hubs_inv4m_cis.csv"))
}

##    module            gene locus_label   kWithin      logFC in_Inv4m
## 1    ME12 Zm00001eb188610        <NA> 4.5743323 -5.9771086    FALSE
## 2    ME12 Zm00001eb188420        <NA> 4.5743323 -6.8786411    FALSE
## 3    ME12 Zm00001eb188430        <NA> 4.5743323 -6.8786411    FALSE
## 4    ME12 Zm00001eb192880        tet2 4.4522112 -6.9457261     TRUE
## 5     ME1 Zm00001eb192530        <NA> 3.5168297  0.5052828     TRUE
## 6     ME1 Zm00001eb195000        <NA> 3.3884404 -0.3487955    FALSE
## 7     ME7 Zm00001eb186400        <NA> 3.3497706  7.7457413    FALSE
## 8     ME1 Zm00001eb193840        <NA> 3.2291615 -0.5468399     TRUE
## 9     ME7 Zm00001eb186390        <NA> 3.0504198  5.2882507    FALSE
## 10    ME4 Zm00001eb190110        <NA> 2.9558591 -0.3443905    FALSE
## 11    ME1 Zm00001eb191080        <NA> 2.9296961  0.4173646     TRUE
## 12    ME1 Zm00001eb191380        <NA> 2.5665566 -0.5453458     TRUE
## 13    ME1 Zm00001eb191990        pig3 2.5517685  1.9529439     TRUE
## 14   ME12 Zm00001eb188460       dof43 2.5093080 -5.1335736    FALSE
## 15    ME4 Zm00001eb195030        <NA> 2.4144039  4.4482408    FALSE
## 16    ME1 Zm00001eb194950        <NA> 2.4079665 -0.4457779    FALSE
## 17    ME6 Zm00001eb192310        <NA> 2.3318359 -0.4493898     TRUE
## 18    ME6 Zm00001eb188040        <NA> 2.2910721 -0.5558307    FALSE
## 19    ME7 Zm00001eb186770      gras80 2.1872245  1.3572539    FALSE
## 20    ME1 Zm00001eb192900        <NA> 2.0397012 -0.8713048     TRUE
## 21    ME7 Zm00001eb186410        <NA> 1.9831002  4.9770347    FALSE
## 22    ME1 Zm00001eb189120        <NA> 1.9485396 -0.5983260    FALSE
## 23    ME1 Zm00001eb190620        <NA> 1.8098399  0.5535087     TRUE
## 24    ME4 Zm00001eb197320        <NA> 1.7985674 -1.1460387    FALSE
## 25    ME7 Zm00001eb186790        adf7 1.7212075  3.6088351    FALSE
## 26   ME12 Zm00001eb196760        <NA> 1.6342953 -8.5057276    FALSE
## 27    ME3 Zm00001eb189030       tola1 1.5654221 -0.4167695    FALSE
## 28    ME1 Zm00001eb188720        <NA> 1.5510188 -0.2309395    FALSE
## 29    ME4 Zm00001eb193790         tu1 1.5143964 -0.4217700     TRUE
## 30    ME2 Zm00001eb196100        abf2 1.5059538  1.7886406    FALSE
## 31    ME5 Zm00001eb197120        <NA> 1.4792710 -0.8027816    FALSE
## 32    ME1 Zm00001eb188530        <NA> 1.4558796  0.3425578    FALSE
## 33    ME1 Zm00001eb191680        <NA> 1.3680893 -0.3345883     TRUE
## 34    ME1 Zm00001eb192340        vps3 1.3631534 -0.6099370     TRUE
## 35    ME1 Zm00001eb192780        <NA> 1.3350413 -0.4478335     TRUE
## 36    ME4 Zm00001eb190150        <NA> 1.3208462 -1.0198821    FALSE
## 37    ME1 Zm00001eb195280        <NA> 1.3097273 -0.2827958    FALSE
## 38    ME2 Zm00001eb193030        <NA> 1.2897970  0.2651506     TRUE
## 39    ME5 Zm00001eb189140        <NA> 1.2787723 -1.1740956    FALSE
## 40    ME1 Zm00001eb196900        <NA> 1.2683280  1.3223904    FALSE
## 41   ME12 Zm00001eb188570      his2b5 1.2564178 -5.3380982    FALSE
## 42    ME5 Zm00001eb190920      nanmt1 1.1891775  0.3457271     TRUE
## 43    ME1 Zm00001eb196500        <NA> 1.1311407 -1.4366628    FALSE
## 44   ME12 Zm00001eb192170        <NA> 1.1196161 -7.4688231     TRUE
## 45    ME4 Zm00001eb192910     lkrsdh1 1.0435040 -1.8177127     TRUE
## 46    ME2 Zm00001eb193530        <NA> 0.9992925 -0.3054514     TRUE
## 47    ME2 Zm00001eb188310        <NA> 0.9896982  0.9907331    FALSE
## 48    ME7 Zm00001eb186610        <NA> 0.9506719  3.9380702    FALSE
## 49    ME3 Zm00001eb193640        <NA> 0.9489488  0.8931164     TRUE
## 50    ME8 Zm00001eb191370       iaa18 0.9317701  0.3571391     TRUE
## 51    ME7 Zm00001eb189850        <NA> 0.9229845 -0.4890672    FALSE
## 52    ME4 Zm00001eb194330        <NA> 0.9205948 -0.6346064     TRUE
## 53    ME4 Zm00001eb193820        <NA> 0.9203769  0.4699150     TRUE
## 54    ME2 Zm00001eb195340        <NA> 0.9182327  1.0326064    FALSE
## 55    ME4 Zm00001eb190090        <NA> 0.9177928  7.1778814    FALSE
## 56    ME4 Zm00001eb189890        <NA> 0.9084722 -1.0289273    FALSE
## 57    ME1 Zm00001eb196360        <NA> 0.8984926  0.2602678    FALSE
## 58    ME1 Zm00001eb190350        <NA> 0.8981436 -1.8752894    FALSE
## 59   ME11 Zm00001eb187910        <NA> 0.8807480 -0.9566888    FALSE
## 60    ME2 Zm00001eb187200        <NA> 0.8506035  0.5562427    FALSE
## 61    ME2 Zm00001eb196060        <NA> 0.8454635 -0.3284765    FALSE
## 62    ME8 Zm00001eb191020     nactf26 0.8078752  0.2409455     TRUE
## 63    ME6 Zm00001eb192850        <NA> 0.7990997 -1.9946674     TRUE
## 64   ME10 Zm00001eb190360        <NA> 0.7498262 -0.9146507    FALSE
## 65    ME7 Zm00001eb186860        <NA> 0.7472795  3.0258799    FALSE
## 66    ME1 Zm00001eb192060        <NA> 0.7242752 -0.3905260     TRUE
## 67    ME5 Zm00001eb192960      ereb14 0.7241475 -0.4713370     TRUE
## 68   ME11 Zm00001eb196630       ufgt4 0.7214190  3.9050799    FALSE
## 69    ME1 Zm00001eb187890         smn 0.6805194  3.9673651    FALSE
## 70    ME7 Zm00001eb186980        <NA> 0.6516490  0.7492750    FALSE
## 71   ME18 Zm00001eb187460        <NA> 0.6486416  2.1857894    FALSE
## 72    ME4 Zm00001eb195270    bnl10.05 0.6468154  0.6908783    FALSE
## 73    ME7 Zm00001eb189210        <NA> 0.6346461 -0.3135679    FALSE
## 74    ME4 Zm00001eb191500        <NA> 0.5984473  0.3208739     TRUE
## 75   ME11 Zm00001eb196090      pmei47 0.5979830  5.3555573    FALSE

Cis Enrichment Analysis

Combined Network

total_cis_all <- sum(module_deg_stats_all$in_cis, na.rm = TRUE)
total_genes_all <- nrow(module_deg_stats_all)

module_cis_enrichment_all <- module_deg_stats_all %>%
  group_by(module) %>%
  summarise(
    n_genes = n(),
    n_in_cis = sum(in_cis, na.rm = TRUE),
    .groups = "drop"
  ) %>%
  rowwise() %>%
  mutate(
    fisher_p = {
      a <- n_in_cis
      b <- n_genes - n_in_cis
      c <- total_cis_all - a
      d <- (total_genes_all - total_cis_all) - b
      
      contingency <- matrix(c(a, b, c, d), nrow = 2)
      fisher.test(contingency)$p.value
    },
    odds_ratio = {
      a <- n_in_cis
      b <- n_genes - n_in_cis
      c <- total_cis_all - a
      d <- (total_genes_all - total_cis_all) - b
      
      (a * d) / (b * c)
    }
  ) %>%
  ungroup() %>%
  mutate(
    fisher_fdr = p.adjust(fisher_p, method = "fdr")
  ) %>%
  arrange(fisher_p)

cat("Cis enrichment analysis (all samples):\n")

## Cis enrichment analysis (all samples):

print(module_cis_enrichment_all)

## # A tibble: 14 × 6
##    module n_genes n_in_cis fisher_p odds_ratio fisher_fdr
##    <chr>    <int>    <int>    <dbl>      <dbl>      <dbl>
##  1 ME1        340      213 7.30e-34     6.63     1.02e-32
##  2 ME2         77        3 1.19e-13     0.0546   8.35e-13
##  3 ME3         60        5 5.71e- 8     0.129    2.67e- 7
##  4 ME4         50        4 6.69e- 7     0.126    2.34e- 6
##  5 ME7         31        6 2.40e- 2     0.366    6.73e- 2
##  6 ME10        29        6 5.10e- 2     0.400    1.19e- 1
##  7 ME14        10        1 9.83e- 2     0.173    1.97e- 1
##  8 ME5         35        9 1.13e- 1     0.532    1.98e- 1
##  9 ME9         29        8 2.47e- 1     0.590    3.84e- 1
## 10 ME13        13        3 3.90e- 1     0.468    5.46e- 1
## 11 ME11        21       10 4.97e- 1     1.45     6.32e- 1
## 12 ME12        20        9 6.44e- 1     1.30     7.51e- 1
## 13 ME6         31       13 7.11e- 1     1.15     7.65e- 1
## 14 ME8         30       11 8.51e- 1     0.910    8.51e- 1

cat("\nModules enriched for cis genes (FDR < 0.05):\n")

## 
## Modules enriched for cis genes (FDR < 0.05):

cis_enriched_all <- module_cis_enrichment_all %>%
  filter(fisher_fdr < 0.05)

if (nrow(cis_enriched_all) > 0) {
  print(cis_enriched_all)
} else {
  cat("  None\n")
}

## # A tibble: 4 × 6
##   module n_genes n_in_cis fisher_p odds_ratio fisher_fdr
##   <chr>    <int>    <int>    <dbl>      <dbl>      <dbl>
## 1 ME1        340      213 7.30e-34     6.63     1.02e-32
## 2 ME2         77        3 1.19e-13     0.0546   8.35e-13
## 3 ME3         60        5 5.71e- 8     0.129    2.67e- 7
## 4 ME4         50        4 6.69e- 7     0.126    2.34e- 6

Genotype-Specific Networks

total_cis_ctrl <- sum(module_deg_stats_ctrl$in_cis, na.rm = TRUE)
total_genes_ctrl <- nrow(module_deg_stats_ctrl)

module_cis_enrichment_ctrl <- module_deg_stats_ctrl %>%
  group_by(module) %>%
  summarise(
    n_genes = n(),
    n_in_cis = sum(in_cis, na.rm = TRUE),
    .groups = "drop"
  ) %>%
  rowwise() %>%
  mutate(
    fisher_p = {
      a <- n_in_cis
      b <- n_genes - n_in_cis
      c <- total_cis_ctrl - a
      d <- (total_genes_ctrl - total_cis_ctrl) - b
      
      contingency <- matrix(c(a, b, c, d), nrow = 2)
      fisher.test(contingency)$p.value
    },
    odds_ratio = {
      a <- n_in_cis
      b <- n_genes - n_in_cis
      c <- total_cis_ctrl - a
      d <- (total_genes_ctrl - total_cis_ctrl) - b
      
      (a * d) / (b * c)
    }
  ) %>%
  ungroup() %>%
  mutate(
    fisher_fdr = p.adjust(fisher_p, method = "fdr")
  ) %>%
  arrange(fisher_p)

total_cis_inv4m <- sum(module_deg_stats_inv4m$in_cis, na.rm = TRUE)
total_genes_inv4m <- nrow(module_deg_stats_inv4m)

module_cis_enrichment_inv4m <- module_deg_stats_inv4m %>%
  group_by(module) %>%
  summarise(
    n_genes = n(),
    n_in_cis = sum(in_cis, na.rm = TRUE),
    .groups = "drop"
  ) %>%
  rowwise() %>%
  mutate(
    fisher_p = {
      a <- n_in_cis
      b <- n_genes - n_in_cis
      c <- total_cis_inv4m - a
      d <- (total_genes_inv4m - total_cis_inv4m) - b
      
      contingency <- matrix(c(a, b, c, d), nrow = 2)
      fisher.test(contingency)$p.value
    },
    odds_ratio = {
      a <- n_in_cis
      b <- n_genes - n_in_cis
      c <- total_cis_inv4m - a
      d <- (total_genes_inv4m - total_cis_inv4m) - b
      
      (a * d) / (b * c)
    }
  ) %>%
  ungroup() %>%
  mutate(
    fisher_fdr = p.adjust(fisher_p, method = "fdr")
  ) %>%
  arrange(fisher_p)

cat("\nCis enrichment analysis (CTRL):\n")

## 
## Cis enrichment analysis (CTRL):

print(module_cis_enrichment_ctrl)

## # A tibble: 19 × 6
##    module n_genes n_in_cis      fisher_p odds_ratio   fisher_fdr
##    <chr>    <int>    <int>         <dbl>      <dbl>        <dbl>
##  1 ME1        116       17 0.00000000183      0.227 0.0000000347
##  2 ME5         56       33 0.00166            2.42  0.0158      
##  3 ME3         82       20 0.00556            0.474 0.0352      
##  4 ME9         36       22 0.00774            2.60  0.0368      
##  5 ME6         40       23 0.0186             2.23  0.0707      
##  6 ME12        27        5 0.0278             0.348 0.0819      
##  7 ME4         63       33 0.0302             1.83  0.0819      
##  8 ME0         27       16 0.0423             2.37  0.100       
##  9 ME8         36       20 0.0525             2.04  0.111       
## 10 ME18        13        2 0.0920             0.282 0.175       
## 11 ME11        28       15 0.116              1.86  0.183       
## 12 ME15        16        3 0.122              0.358 0.183       
## 13 ME10        30       16 0.125              1.85  0.183       
## 14 ME14        17        4 0.219              0.479 0.297       
## 15 ME7         39       12 0.317              0.689 0.401       
## 16 ME16        15        7 0.596              1.39  0.708       
## 17 ME2        102       41 0.745              1.07  0.833       
## 18 ME17        15        5 0.792              0.785 0.836       
## 19 ME13        18        7 1                  1.00  1

cat("\nModules enriched for cis genes (CTRL, FDR < 0.05):\n")

## 
## Modules enriched for cis genes (CTRL, FDR < 0.05):

cis_enriched_ctrl <- module_cis_enrichment_ctrl %>%
  filter(fisher_fdr < 0.05)

if (nrow(cis_enriched_ctrl) > 0) {
  print(cis_enriched_ctrl)
} else {
  cat("  None\n")
}

## # A tibble: 4 × 6
##   module n_genes n_in_cis      fisher_p odds_ratio   fisher_fdr
##   <chr>    <int>    <int>         <dbl>      <dbl>        <dbl>
## 1 ME1        116       17 0.00000000183      0.227 0.0000000347
## 2 ME5         56       33 0.00166            2.42  0.0158      
## 3 ME3         82       20 0.00556            0.474 0.0352      
## 4 ME9         36       22 0.00774            2.60  0.0368

cat("\nCis enrichment analysis (Inv4m):\n")

## 
## Cis enrichment analysis (Inv4m):

print(module_cis_enrichment_inv4m)

## # A tibble: 19 × 6
##    module n_genes n_in_cis  fisher_p odds_ratio fisher_fdr
##    <chr>    <int>    <int>     <dbl>      <dbl>      <dbl>
##  1 ME12        16       14 0.0000730     11.5      0.00139
##  2 ME3         45        7 0.000776       0.274    0.00738
##  3 ME1         62       35 0.00405        2.18     0.0257 
##  4 ME6         38        7 0.00959        0.341    0.0456 
##  5 ME14        14        1 0.0124         0.118    0.0473 
##  6 ME4         44       25 0.0161         2.17     0.0508 
##  7 ME17        12        1 0.0346         0.141    0.0940 
##  8 ME7         33       18 0.0680         1.95     0.161  
##  9 ME9         21        5 0.178          0.485    0.346  
## 10 ME13        15        3 0.182          0.388    0.346  
## 11 ME10        19        5 0.343          0.556    0.564  
## 12 ME2         50       16 0.369          0.728    0.564  
## 13 ME16        12        3 0.386          0.521    0.564  
## 14 ME8         26       12 0.540          1.37     0.732  
## 15 ME15        13        6 0.579          1.36     0.734  
## 16 ME0        288      115 0.647          1.08     0.769  
## 17 ME18        11        5 0.758          1.32     0.847  
## 18 ME11        17        7 0.808          1.11     0.852  
## 19 ME5         40       16 0.869          1.05     0.869

cat("\nModules enriched for cis genes (Inv4m, FDR < 0.05):\n")

## 
## Modules enriched for cis genes (Inv4m, FDR < 0.05):

cis_enriched_inv4m <- module_cis_enrichment_inv4m %>%
  filter(fisher_fdr < 0.05)

if (nrow(cis_enriched_inv4m) > 0) {
  print(cis_enriched_inv4m)
} else {
  cat("  None\n")
}

## # A tibble: 5 × 6
##   module n_genes n_in_cis  fisher_p odds_ratio fisher_fdr
##   <chr>    <int>    <int>     <dbl>      <dbl>      <dbl>
## 1 ME12        16       14 0.0000730     11.5      0.00139
## 2 ME3         45        7 0.000776       0.274    0.00738
## 3 ME1         62       35 0.00405        2.18     0.0257 
## 4 ME6         38        7 0.00959        0.341    0.0456 
## 5 ME14        14        1 0.0124         0.118    0.0473

Match Trans Network Edges to WGCNA Modules

cat("Loading trans co-expression network...\n")

## Loading trans co-expression network...

trans_edges <- read_csv(file.path(paths$intermediate, "inv4m_bipartite_edges.csv"))

cat("Trans network edges loaded:", nrow(trans_edges), "\n")

## Trans network edges loaded: 3766

cat("Columns:", paste(colnames(trans_edges), collapse = ", "), "\n")

## Columns: from, to, r

module_lookup_all <- mod_all %>%
  select(gene, module)

trans_edges$in_wgcna_all <- mapply(
  function(f, t) {
    exists_forward <- any(edge_list_all$from == f & edge_list_all$to == t)
    exists_reverse <- any(edge_list_all$from == t & edge_list_all$to == f)
    
    exists_forward | exists_reverse
  },
  trans_edges$from,
  trans_edges$to
)

trans_edges$wgcna_weight_all <- mapply(
  function(f, t) {
    forward_match <- edge_list_all %>%
      filter(from == f & to == t) %>%
      pull(weight)
    
    reverse_match <- edge_list_all %>%
      filter(from == t & to == f) %>%
      pull(weight)
    
    if (length(forward_match) > 0) {
      return(forward_match[1])
    } else if (length(reverse_match) > 0) {
      return(reverse_match[1])
    } else {
      return(NA_real_)
    }
  },
  trans_edges$from,
  trans_edges$to
)

trans_edges_all <- trans_edges %>%
  left_join(
    module_lookup_all,
    by = c("to" = "gene")
  ) %>%
  rename(module_all = module)

cat("\nTrans edges in WGCNA network (all samples):", 
    sum(trans_edges$in_wgcna_all), "out of", nrow(trans_edges), "\n")

## 
## Trans edges in WGCNA network (all samples): 3766 out of 3766

trans_summary_all <- trans_edges_all %>%
  filter(!is.na(module_all)) %>%
  group_by(module_all) %>%
  summarise(
    n_trans_edges = n(),
    n_in_wgcna = sum(in_wgcna_all, na.rm = TRUE),
    pct_in_wgcna = 100 * n_in_wgcna / n_trans_edges,
    n_unique_trans_genes = n_distinct(to),
    n_unique_cis_regulators = n_distinct(from),
    mean_abs_r = mean(abs(r)),
    mean_wgcna_weight = mean(wgcna_weight_all, na.rm = TRUE),
    .groups = "drop"
  ) %>%
  arrange(desc(n_trans_edges))

cat("\nTrans edges per WGCNA module (all samples):\n")

## 
## Trans edges per WGCNA module (all samples):

print(trans_summary_all)

## # A tibble: 5 × 8
##   module_all n_trans_edges n_in_wgcna pct_in_wgcna n_unique_trans_genes
##   <chr>              <int>      <int>        <dbl>                <int>
## 1 ME1                 3473       3473          100                   33
## 2 ME8                   89         89          100                    1
## 3 ME14                  73         73          100                    1
## 4 ME12                  68         68          100                    1
## 5 ME3                   63         63          100                    1
## # ℹ 3 more variables: n_unique_cis_regulators <int>, mean_abs_r <dbl>,
## #   mean_wgcna_weight <dbl>

write_csv(trans_edges_all, file.path(paths$intermediate, "trans_modules_all.csv"))
write_csv(trans_summary_all, file.path(paths$intermediate, "trans_summary_all.csv"))

cat("\nExported:\n")

## 
## Exported:

cat("  trans_modules_all.csv\n")

##   trans_modules_all.csv

cat("  trans_summary_all.csv\n")

##   trans_summary_all.csv

Adjacency Matrix Visualizations

Combined Network (All Samples)

cat("Calculating TOM (all samples)...\n")

## Calculating TOM (all samples)...

TOM_all <- TOMsimilarity(adj_all)

## ..connectivity..
## ..matrix multiplication (system BLAS)..
## ..normalization..
## ..done.

rownames(TOM_all) <- rownames(adj_all)
colnames(TOM_all) <- colnames(adj_all)

cat("Using WGCNA dendrogram ordering...\n")

## Using WGCNA dendrogram ordering...

all_genes <- colnames(TOM_all)
genes_with_modules <- names(net_all$colors)[net_all$colors != 0]

TOM_all_filtered <- TOM_all[genes_with_modules, genes_with_modules]

if (!is.null(net_all$dendrograms) && length(net_all$dendrograms) > 0) {
  full_order <- net_all$dendrograms[[1]]$order
  all_genes_ordered <- all_genes[full_order]
  genes_in_modules_ordered <- all_genes_ordered[all_genes_ordered %in% genes_with_modules]
  gene_order_all <- match(genes_in_modules_ordered, colnames(TOM_all_filtered))
} else {
  gene_order_all <- order(net_all$colors[genes_with_modules], genes_with_modules)
}

TOM_all_ordered <- TOM_all_filtered[gene_order_all, gene_order_all]
TOM_all_ordered_rev <- TOM_all_ordered[nrow(TOM_all_ordered):1, ]

gene_names_ordered_all <- colnames(TOM_all_filtered)[gene_order_all]
gene_names_ordered_all_rev <- rev(gene_names_ordered_all)

gene_module_annotation_all <- mod_all %>%
  filter(gene %in% gene_names_ordered_all_rev) %>%
  select(gene, module) %>%
  mutate(gene = factor(gene, levels = gene_names_ordered_all_rev)) %>%
  arrange(gene) %>%
  column_to_rownames("gene")

pheatmap(
  TOM_all_ordered_rev,
  annotation_row = gene_module_annotation_all,
  annotation_col = gene_module_annotation_all,
  cluster_rows = FALSE,
  cluster_cols = FALSE,
  show_rownames = FALSE,
  show_colnames = FALSE,
  main = paste0("WGCNA TOM Similarity - All Samples\n(power = ", 
                pow_all, ")"),
  breaks = seq(0, 1, length.out = 101),
  color = colorRampPalette(c("white", "red"))(100)
)

write_csv(
  as.data.frame(TOM_all_ordered_rev) %>% rownames_to_column("gene"),
  file.path(paths$intermediate, "tom_all.csv")
)

cat("\nWGCNA TOM matrix exported (all samples)\n")

## 
## WGCNA TOM matrix exported (all samples)

Genotype-Specific Split TOM Matrices

cat("Calculating TOM (genotype-specific)...\n")

## Calculating TOM (genotype-specific)...

TOM_ctrl <- TOMsimilarity(adj_ctrl)

## ..connectivity..
## ..matrix multiplication (system BLAS)..
## ..normalization..
## ..done.

rownames(TOM_ctrl) <- rownames(adj_ctrl)
colnames(TOM_ctrl) <- colnames(adj_ctrl)

TOM_inv4m <- TOMsimilarity(adj_inv4m)

## ..connectivity..
## ..matrix multiplication (system BLAS)..
## ..normalization..
## ..done.

rownames(TOM_inv4m) <- rownames(adj_inv4m)
colnames(TOM_inv4m) <- colnames(adj_inv4m)

cat("TOM matrices calculated\n")

## TOM matrices calculated

all_genes_inv4m <- colnames(TOM_inv4m)
genes_with_modules_inv4m <- names(net_inv4m$colors)[net_inv4m$colors != 0]

TOM_inv4m_filtered <- TOM_inv4m[genes_with_modules_inv4m, genes_with_modules_inv4m]
TOM_ctrl_filtered <- TOM_ctrl[genes_with_modules_inv4m, genes_with_modules_inv4m]

cat("Using Inv4m WGCNA dendrogram ordering...\n")

## Using Inv4m WGCNA dendrogram ordering...

if (!is.null(net_inv4m$dendrograms) && length(net_inv4m$dendrograms) > 0) {
  full_order_inv4m <- net_inv4m$dendrograms[[1]]$order
  all_genes_ordered_inv4m <- all_genes_inv4m[full_order_inv4m]
  genes_in_modules_ordered_inv4m <- all_genes_ordered_inv4m[all_genes_ordered_inv4m %in% genes_with_modules_inv4m]
  gene_order_inv4m <- match(genes_in_modules_ordered_inv4m, colnames(TOM_inv4m_filtered))
} else {
  gene_order_inv4m <- order(net_inv4m$colors[genes_with_modules_inv4m], genes_with_modules_inv4m)
}

TOM_ctrl_ordered_inv4m <- TOM_ctrl_filtered[gene_order_inv4m, gene_order_inv4m]
TOM_inv4m_ordered_inv4m <- TOM_inv4m_filtered[gene_order_inv4m, gene_order_inv4m]

combined_tom_inv4m <- TOM_ctrl_ordered_inv4m
combined_tom_inv4m[lower.tri(combined_tom_inv4m)] <- 
  TOM_inv4m_ordered_inv4m[lower.tri(TOM_inv4m_ordered_inv4m)]

combined_tom_inv4m_t <- t(combined_tom_inv4m)
combined_tom_inv4m_rev <- combined_tom_inv4m_t[nrow(combined_tom_inv4m_t):1, ]

gene_names_ordered_inv4m <- colnames(TOM_inv4m_filtered)[gene_order_inv4m]
gene_names_ordered_inv4m_rev <- rev(gene_names_ordered_inv4m)

gene_module_annotation_inv4m <- mod_inv4m %>%
  filter(gene %in% gene_names_ordered_inv4m_rev) %>%
  select(gene, module) %>%
  mutate(gene = factor(gene, levels = gene_names_ordered_inv4m_rev)) %>%
  arrange(gene) %>%
  column_to_rownames("gene")

pheatmap(
  combined_tom_inv4m_rev,
  annotation_row = gene_module_annotation_inv4m,
  annotation_col = gene_module_annotation_inv4m,
  cluster_rows = FALSE,
  cluster_cols = FALSE,
  show_rownames = FALSE,
  show_colnames = FALSE,
  main = paste0("WGCNA TOM Similarity (Upper: CTRL [pow=", pow_ctrl, 
                "], Lower: Inv4m [pow=", pow_inv4m, "])\nClustered by Inv4m"),
  breaks = seq(0, 1, length.out = 101),
  color = colorRampPalette(c("white", "red"))(100)
)

write_csv(
  as.data.frame(TOM_ctrl_ordered_inv4m) %>% rownames_to_column("gene"),
  file.path(paths$intermediate, "tom_ctrl_inv4m_order.csv")
)

write_csv(
  as.data.frame(TOM_inv4m_ordered_inv4m) %>% rownames_to_column("gene"),
  file.path(paths$intermediate, "tom_inv4m_inv4m_order.csv")
)

all_genes_ctrl <- colnames(TOM_ctrl)
genes_with_modules_ctrl <- names(net_ctrl$colors)[net_ctrl$colors != 0]

TOM_ctrl_filtered_ctrl <- TOM_ctrl[genes_with_modules_ctrl, genes_with_modules_ctrl]
TOM_inv4m_filtered_ctrl <- TOM_inv4m[genes_with_modules_ctrl, genes_with_modules_ctrl]

cat("Using CTRL WGCNA dendrogram ordering...\n")

## Using CTRL WGCNA dendrogram ordering...

if (!is.null(net_ctrl$dendrograms) && length(net_ctrl$dendrograms) > 0) {
  full_order_ctrl <- net_ctrl$dendrograms[[1]]$order
  all_genes_ordered_ctrl <- all_genes_ctrl[full_order_ctrl]
  genes_in_modules_ordered_ctrl <- all_genes_ordered_ctrl[all_genes_ordered_ctrl %in% genes_with_modules_ctrl]
  gene_order_ctrl <- match(genes_in_modules_ordered_ctrl, colnames(TOM_ctrl_filtered_ctrl))
} else {
  gene_order_ctrl <- order(net_ctrl$colors[genes_with_modules_ctrl], genes_with_modules_ctrl)
}

TOM_ctrl_ordered_ctrl <- TOM_ctrl_filtered_ctrl[gene_order_ctrl, gene_order_ctrl]
TOM_inv4m_ordered_ctrl <- TOM_inv4m_filtered_ctrl[gene_order_ctrl, gene_order_ctrl]

combined_tom_ctrl <- TOM_inv4m_ordered_ctrl
combined_tom_ctrl[lower.tri(combined_tom_ctrl)] <- 
  TOM_ctrl_ordered_ctrl[lower.tri(TOM_ctrl_ordered_ctrl)]

combined_tom_ctrl_t <- t(combined_tom_ctrl)
combined_tom_ctrl_rev <- combined_tom_ctrl_t[nrow(combined_tom_ctrl_t):1, ]

gene_names_ordered_ctrl <- colnames(TOM_ctrl_filtered_ctrl)[gene_order_ctrl]
gene_names_ordered_ctrl_rev <- rev(gene_names_ordered_ctrl)

gene_module_annotation_ctrl <- mod_ctrl %>%
  filter(gene %in% gene_names_ordered_ctrl_rev) %>%
  select(gene, module) %>%
  mutate(gene = factor(gene, levels = gene_names_ordered_ctrl_rev)) %>%
  arrange(gene) %>%
  column_to_rownames("gene")

pheatmap(
  combined_tom_ctrl_rev,
  annotation_row = gene_module_annotation_ctrl,
  annotation_col = gene_module_annotation_ctrl,
  cluster_rows = FALSE,
  cluster_cols = FALSE,
  show_rownames = FALSE,
  show_colnames = FALSE,
  main = paste0("WGCNA TOM Similarity (Upper: Inv4m [pow=", pow_inv4m, 
                "], Lower: CTRL [pow=", pow_ctrl, "])\nClustered by CTRL"),
  breaks = seq(0, 1, length.out = 101),
  color = colorRampPalette(c("white", "red"))(100)
)

write_csv(
  as.data.frame(TOM_ctrl_ordered_ctrl) %>% rownames_to_column("gene"),
  file.path(paths$intermediate, "tom_ctrl_ctrl_order.csv")
)

write_csv(
  as.data.frame(TOM_inv4m_ordered_ctrl) %>% rownames_to_column("gene"),
  file.path(paths$intermediate, "tom_inv4m_ctrl_order.csv")
)

cat("\nGenotype-specific TOM matrices exported\n")

## 
## Genotype-specific TOM matrices exported